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Abstract. The engineering and control of devices at the quantum-mechanical level — such as 
those consisting of small numbers of atoms and photons — is a delicate business. The fundamental 
uncertainty that is inherently present at this scale manifests itself in the unavoidable presence of noise, 
making this a novel field of application for stochastic estimation and control theory. In this expository 
paper we demonstrate estimation and feedback control of quantum mechanical systems in what is 
essentially a noncommutative version of the binomial model that is popular in mathematical finance. 
The model is extremely rich and allows a full development of the theory, while remaining completely 
within the setting of finite-dimensional Hilbert spaces (thus avoiding the technical complications 
of the continuous theory). We introduce discretized models of an atom in interaction with the 
electromagnetic field, obtain filtering equations for photon counting and homodyne detection, and 
solve a stochastic control problem using dynamic programming and Lyapunov function methods. 

Key words, discrete quantum filtering, quantum feedback control, quantum probability, con- 
ditional expectation, dynamic programming, stochastic Lyapunov functions. 

AMS subject classifications. 93E11, 93E15, 93E20, 81P15, 81S25, 34F05. 

Dedicated to Slava Belavkin in the year of his 60th birthday. 

1. Introduction. Control theory, and in particular feedback control, is an im- 
portant aspect of modern engineering. It provides a set of tools for the design of 
technology with reliable performance and has been applied with overwhelming suc- 
cess in the design of many of the devices that we use on a daily basis. In this article 
we will explore the following question: rather than controlling, say, a jet engine, can 
we use feedback to control an object as small as a single atom? 

Though it is not directly visible in our everyday lives, the technology to manipu- 
late matter at the level of single atoms and photons is well in place and is progressing 
at a fast rate. Nobel prize winning technologies such as laser cooling and trapping 
of atoms, which were state-of-the-art only a decade ago, are now routine procedures 
which are being applied in physics laboratories around the world. The generation of 
highly coherent light and the detection of single photons has been refined since the 
development of the laser. Given the availability of this unprecedented level of control 
over the microscopic world, the question of control design seems a timely one — one 
could wonder whether we can "close the loop" and make the atoms do useful work 
for us. This is the domain of quantum feedback control. 

"Quantum," of course, refers to the theory of quantum mechanics, which comes 
in necessarily at this level of description. The feedback which we will use to control 
a quantum system is a function of observations obtained from that system, and ob- 
servations in quantum mechanics are inherently random. This makes the theory both 
challenging and interesting from the point of view of fundamental physics and control 
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Fig. 1.1. A laboratory experiment in quantum feedback control, implementing the setup reported 
in |39) . The metal boxes in the foreground are diode lasers, whose output light is manipulated using 
a large number of optical elements. In the background on the left is a magnetically shielded chamber 
containing magnetic coils and a cell with laser cooled Cesium atoms. A schematic of the experiment 
is shown on the right. (Photograph courtesy of John Stockton.) 

theory. As we will demonstrate, the theory of quantum feedback control resembles 
closely, and follows directly in many parts, the classical theory of stochastic control 
(here and below, classical means non-quantum mechanical). 

Several promising applications of quantum feedback control have been proposed 
and are now starting to be investigated in a laboratory setting. One class of such 
applications falls under the general heading of precision metrology: can we utilize the 
sensitivity of small numbers of atoms to external perturbations to design ultrafast, 
precise sensors, when the desired accuracy is on the order of the intrinsic quantum 
uncertainty of the sensor? Concrete applications include e.g. precision magnetometry 
[38l l77] , which was recently investigated in a laboratory experiment [39] (see Figure 
ll.ip . and atomic clocks 2J. A second class of examples involves the design of optical 
communication systems where each bit is encoded in a pulse containing only a small 
number of photons [551 [37]; see [3] for a related experimental demonstration. As a 
third application, we mention the use of feedback for cooling of atoms or nanome- 
chanical devices [311 \M\ [IS] , see [23] for a laboratory demonstration (the availability 
of efficient cooling techniques is important in experimental physics). 

As both the theory and applications are still in their infancy, it is far from clear at 
this time what the ultimate impact of quantum feedback control will be; it is clearly 
unlikely that future generations of washing machines or airplanes will be powered by 
quantum control, but other application areas, such as the ones mentioned above, might 
well benefit from the availability of such technology. At the very least, however, the 
questions that arise in the subject prompt us to explore the fundamental limits on the 
control and engineering of physical systems, and a fruitful interplay between control 
theory and experimental physics enables the exploration of these limits using real-life 
experiments such as the one displayed in Figure 11.11 In this article we demonstrate 
the underlying theory using a simple but rich toy model; we hope that the reader is 
sufficiently intrigued to come along for the ride! 

The study of quantum feedback control was pioneered by V. P. Belavkin in a 1983 
paper [S], long before the experimental realization of such a setup was realistically 
feasible. Belavkin realized that due to the unavoidable presence of uncertainty in 
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quantum measurements, the theory of filtering — the extraction of information from 
noisy observations — plays a fundamental role in quantum feedback control, exactly 
as in the theory of stochastic control with partial observations [T^. In [Hj feedback 
control was explored in discrete time using the operational formalism of Davies [27j , 
which is the precursor of quantum filtering theory. Belavkin developed further the 
theory of quantum filtering in a series of articles [SI [TOl [12] , and in [10] the use of 
dynamic programming for feedback control in continuous time was sketched. We refer 
to [20j for a modern introduction to quantum filtering theory. 

Independently from Belavkin's work, the equations of quantum nonlinear filtering 
were developed in the physics literature by Carmichael [24j . based on the work of 
Davies [27], under the name of "quantum trajectory theory". Though the connection 
to filtering theory was not realized until much later, Carmichael's work nonetheless 
spawned an investigation on the use of feedback control in the physics literature (see 
e.g. [3D] and the references therein). It is here that most of the practical applications, 
such as the ones cited above, were proposed and developed. 

Like its classical counterpart, the mathematical theory of quantum filtering and 
stochastic control in continuous time can be quite technical, requiring a heavy analytic 
machinery to obtain the necessary results in a rigorous way. On the other hand, we 
believe that the methods in the field can be understood perfectly well without being 
obscured by analytic complications. What we are going to do in this paper is to give a 
complete development of the basic theory of quantum filtering and feedback control in 
a simple toy model, which requires essentially no analytic machinery. We only assume 
some background in linear algebra and familiarity with elementary probability theory 
with martingales, roughly at the level of the inspiring textbook by D. Williams [82] . 

The model we will investigate likely has a familiar ring to those readers who are 
familiar with mathematical finance. It is in fact little more than a noncommutative 
version of the binomial model which has enjoyed wide popularity since its introduction 
by Cox, Ross and Rubinstein ^26j. The model is widely used in the teaching of 
mathematical finance [75] and even on Wall Street as a flexible computational tool 
for the pricing of financial derivatives. The model has three key features, all of 
which are shared by its noncommutative counterpart, which make its investigation 
particularly appealing: (i) on the one hand the model is as simple as it can possibly 
be, consisting of a finite number of two-state random variables (coin flips) — hence 
analytical complications are avoided; (ii) on the other hand, many of the mathematical 
tools that are indispensable in the continuous context can be developed within the 
framework of this model and admit simple proofs: in particular, change of measure 
techniques, the martingale representation theorem, discrete stochastic integrals and 
even a (trivial) stochastic calculus can all be demonstrated in this context; and (iii) 
the model is in fact an approximation to the usual continuous time models in the 
sense that these are recovered when the limit as the time step tends to zero is taken 
in a suitable sense. For these reasons, we believe that such models are ideally suited 
as a first introduction to quantum filtering and feedback control. 

The noncommutative model has an independent history in the quantum proba- 
bility literature. In his book [69] P.-A. Meyer uses the building blocks of the model 
we use in this article to demonstrate the basic features of the Fock space, which un- 
derlies the continuous quantum stochastic calculus of Hudson-Parthasarathy [50| , in 
a discrete setting (he calls this model, which he attributes to Journe [69l page 18], 
the "toy Fock space"). Repeated interaction models of a similar type were consid- 
ered by Kiimmerer |62| . Various authors have considered the continuous limit of such 
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models [5U |3 [B] |35] and have demonstrated convergence of the associated discretized 
evolution equations to (Hudson-Parthasarathy) quantum stochastic differential equa- 
tions which are routinely used for the description of realistic physical systems |36| . 
The description of measurements (in essence, filtering) in this discrete framework was 
deduced from basic physical arguments in Brun [22] and Gough [44] , 

In this paper we will not take the simplest and most straightforward route in 
treating these discrete models; similar equations can often be deduced by more ele- 
mentary means |22i 144] . Rather, we will exploit the richness of the model to mimic 
as closely as possible the development of quantum filtering and feedback control in 
the continuous setting, where simple arguments based on linear algebra are useless. 
Indeed, the introduction of a "trivial" stochastic calculus, martingale representation, 
etc., may seem overkill in the context of this model, but these become powerful and 
indispensable tools in the treatment of the more general (and physically relevant) 
framework of continuous time evolution and observations. We hope to provide in 
this way a (relatively) gentle introduction to the full technical theory, where all the 
necessary ingredients make their appearance in a more accessible setting. 

The remainder of the paper is organized as follows. In fj2]we introduce the basic 
notions of quantum probability, and develop in detail the quantum binomial model 
which we will be studying throughout the article. f}3] introduces conditional expec- 
tations in the quantum context and sets up the filtering problem. In Sj4] we develop 
a discrete stochastic calculus and use it to express the model of in the form of a 
difference equation. fj5] treats filtering using martingale methods, while <j6]uses change 
of measure techniques. In Sj7] we show how to incorporate feedback into the model 
of Sj2j In ^J8] we develop optimal feedback controls using dynamic programming tech- 
niques, and fj9] briefly explores an alternative control design method using Lyapunov 
functions. We conclude in fJTUlwith a set of references for further reading. 

Notation. We use the following conventions. In general, Hilbert spaces are 
denoted by the symbol H and are endowed with the inner product (•,•). We will 
assume all Hilbert spaces are complex and finite-dimensional. Recall that the adjoint 
of a linear operator X : H — > H is defined as the unique operator X* that satisfies 
{X*x, y) = (x, Xy) Vx, y G H, and that X is called self-adjoint if X = X* . We denote 
by I the identity operator. The commutator between two linear operators is denoted 
as [X, Y] — XY — YX. Calligraphic letters will be used for various purposes, such as 
classical c-algebras (e.g. y) or functionals on an algebra (e.g. C{X)). Control sets and 
related quantities are denoted by gothic type (e.g. il). *-algebras (to be defined below) 
will be denoted by script type and states on such an algebra are often denoted 
by blackboard type (e.g. P), though p and </) will also be used. The generic classi- 
cal probability space will be denoted as (ri,JF, P), Ep denotes the expectation with 
respect to the measure P, and £°°{!F) denotes the space of JF-measurable (bounded, 
but trivially so when f2 is a finite set) random variables. 

As is unavoidable in an article of this length, there is much notation that is in- 
troduced throughout the article and that is constantly reused. To help the reader 
to keep track of the various quantities, we have provided below a list of commonly 
used objects, a brief description, and the page number on which the object is defined. 



A{l),A*{l),A{l),t{l) Discrete noises 
J(/x) Control cost function 

ii.2,3 Interaction operators 

M{1) Single time step interaction 
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M^^{l,u) Controlled single time step interaction 51 

j\//°,+^-^± CocfBcicnts in difference equation for [/(/) 26 

Ml Single time step interaction, time reversed 16 

Mi{u) Controlled single time step interaction, time reversed 50 

U{1) Interaction unitary 16 

U{l,u) Interaction unitary with open loop control 49 

U^{1) Interaction unitary with feedback control 51 

V{1) Change of state operator for reference measure 43 

Vi{p) Value function of dynamic programming 59 

Wi{fi,p) Cost-to-go 58 

Y{1) Observation process 17 

F'^(Z) Observations process for control strategy /x 53 

Z{1) Field process with Y{1) = U (l)* Z {l)U (l) 17 

AL{1) 24 

$ Vacuum vector in 10 

L *-isomorphism between observables and random variables 8 

Length of each time slice 11 

P Standard state on the atom-field algebra ®Wk 15 

J{X^ jT'-coefRcient of the filtering equations 34 

C{X) Discrete Lindblad generator 27 

T(X) T-coefficient of the filtering equations 47 

yi Classical observation filtration 18 

^ Set of admissible controls 51 

^5 Set of admissible separated controls 56 

it Control set 49 

Ui Control process /,(AZ(1), . . . , AZ(Z - 1)) 52 

U; Control signal /i(Ay(l),...,Ar(/-l)) 51 

3§i Filtration in ^ 24 

% Filtration of ^ 17 

^ Two-level system algebra 9 

Ml Field algebra of time slice / 12 

Wu Field algebra 12 

'3^1 Observation filtration 18 

IJ,,H*,p, Control strategies 51 

(j) Vacuum state on M 13 

7ri{X) Conditional state (filtered estimate of X) 23 

f±)fz Pauli matrices in ^ 10 

(Ji{X) Unnormalized conditional state 44 

Y{1) Innovations process 31 

/;(•••) Feedback function 51 

gi{p) Separated feedback function 56 

ji{X) Time evolution of X 16 

Time evolution with open loop control 49 

jl'iX) Time evolution with feedback control 51 

k Number of time slices 1 1 

uj* Classical feedback signal for control strategy fj, 54 

Classical observations process for control strategy /x 54 

yi Classical observation process 18 
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2. The quantum binomial model. In this section we introduce the basic 
model that we will be dealing with throughout: a discretized approximation of the 
interaction between an atom and the electromagnetic field. First, however, we need 
to demonstrate how probability theory fits in the framework of quantum mechanics. 

2.1. Random variables in quantum mechanics. The basic setting of quan- 
tum mechanics, as one would find it in most textbooks, is something like this. We 
start with a Hilbert space H and fix some "state vector" ip G H. An "observable", 
the physicist's word for random variable, is described by a self-adjoint operator X 
on H, and the expectation of X is given by {ip,Xip). The set of values that X can 
take in a single measurement is its set of eigenvalues, and the probability of observing 
the eigenvalue Xi is given by {tp, Pitp) where Pi is the projection operator onto the 
eigenspace corresponding to A^. This is quite unlike the sort of description we are 
used to from classical probability theory — or is it? 

In fact, the two theories are not as dissimilar as they may seem, and it is fruitful 
to formalize this idea (we will do this in the next section). The key result that we need 
is the following elementary fact from linear algebra. This is just the statement that a 
set of commuting normal ([X, X*] = 0) matrices can be simultaneously diagonalized, 
see e.g. gSJ sec. 2.5] or [45l sec. 84]. 

Theorem 2.1 (Spectral theorem). Let H be an n- dimensional Hilbert space, 
n < oo. Let ^S' be a set of linear transformations from H — > H that is closed under the 
adjoint {i.e. if C £^ then also C* € '^) and such that all the elements of^ commute 
[i.e. [C, D] — "iCjD G '^). Then there exists an orthonormal basis of H such that 
every C £^ is represented by a diagonal matrix with respect to this basis. 

Let us demonstrate how this works in the simplest case. Let dim H = 2, fix some 
■0 g H and let X = X* be a self-adjoint operator on H. The set 'rf = {X} satisfies 
the conditions of the spectral theorem, so we can find an orthonormal basis in H such 
that we can express X and ip in this basis as 

We can now interpret X as a random variable on some probability space. Introduce 

= {1, 2}, the map x : ^ R, x{i) = Xi, and the probability measure P({j}) = 
[i — 1,2). Evidently (V',Xt/;) — Ep(a;), i.e. the random variable x has the same 
expectation under P as we obtain from X and ij} using the usual quantum mechanical 
formalism. We can also easily calculate P(x = Xi) = IV'iP, which is again consistent 
with the rules of quantum mechanics. As the spectral theorem works for more general 
sets we can follow the same procedure to represent a set of commuting observables 
on a classical probability space and to calculate the joint statistics. 

Up to this point we have done nothing particularly exciting: all we have done is to 
express the quantum "rules" listed in the first paragraph of this section in a convenient 
basis, and we have attached an interpretation in terms of classical probability theory. 
Conceptually, however, this is an important point: commuting observables are random 
variables on a classical probability space, and should be thought of in that way. 
Formalizing the way we can pass back and forth between the two descriptions will 
allow us to utilize classical probabilistic techniques in the quantum mechanical setting. 

What the spectral theorem does not allow us to do is to simultaneously interpret 
two noncommuting observables as random variables on a single probability space. 
This is not a shortcoming of the theory, but has an important physical meaning. Ob- 
servables that do not commute cannot be observed in the same realization: their joint 
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probability distribution is a meaningless quantity (hence they are called incompati- 
ble). Strange as this may seem, this idea is a cornerstone of quantum theory and all 
the empirical evidence supports the conclusion that this is how nature works. We will 
accept as a physical principle that in a single realization we can choose to measure 
at most a commiiting set of observables. However, we will see later that we can still 
estimate the statistics of observables which we did not choose to measure, even if these 
do not commute with each other, provided they commute with the measurements we 
did choose to perform. These ideas will be clarified in due cotirse. 

2.2. Quantum probability spaces. In this section we are going to introduce 
and motivate the notion of a (finite-dimensional) quantum probability space. The 
definitions may not seem completely natural at first sight, but their use will be amply 
illustrated in the remainder of the article. 

In the example of the previous section, we saw how to construct a set of quan- 
tum observables as maps on the sample space fl. To turn this into a probability 
space we need to add a u-algebra and a probability measure. As we are considering 
a particular set of random variables Xi , . . . , it is natural to choose the cr-algebra 
<t{Xi, . . . , Xfc), i.e. the smallest cr-algebra with respect to which Xi, . . . ,Xk are mea- 
surable (the (T-algebra generated by Xi, . . . ,Xk), and the quantum state induces a 
probability measure on this a-algebra. We would like to be able to specify the cr- 
algebra directly, however, without having to select some "preferred" set of random 
variables Xi, . . . ,Xk the values taken by these random variables are irrelevant, as 
we are only interested in the generated cr-algebra. 

The definition we will give is motivated by the following observation. On a clas- 
sical sample space ft, consider the set £°°{!F) of all random variables measurable with 
respect to the cr-algebra J^. If one is given such a set, then can be reconstructed 
as .F = a{£°°{!F)}. Conversely, however, we can easily construct if we are 

given J^. Hence these two structm-es carry the same information; by considering all 
measurable random variables at the same time we are no longer giving preference to 
a particular set Xi, . . . ,X}~. In quantum probability, it will be more convenient to 
characterize a cr-algebra by its £°°-space — the quantum analog of this object is a set 
of observables jz/, which can be used directly as input for the spectral theorem. To 
complete the picture we need a way to characterize sets of observables which are 
mapped by the spectral theorem to a space of the form £°°{!F). 

Definition 2.2. A *-algebra is a set of linear transformations H ^ H such 
that I, aA+l3B, AB, A* G £/ for any A,Bg , a,f3 G C is called commutative if 
AB = BA for any A, B G . A linear map p : ^ C which is positive p{A*A) > 
WA G and normalized p{I) = 1 is called a state on si . 

Definition 2.3. A * -isomorphism between a commutative * -algebra and a set 
of functions A on some space fl is a linear bijection l : s/ A such that ^iA*) = 
l{A)* {L{A)*{i) is the com.plex conjugate of i.(A){i) Vi € ft) and l{AB) = L{A)i{B) 
{{L{A)L{B)){i) = L{A){i)L{B){i) Vi e fl) for every A, B € s/ . 

The reason we want *-isomorphisms is to ensure that we can manipulate observ- 
ables and classical random variables in the same way. That is, if Xi, X2 are commuting 
self-adjoint operators that correspond to the random variables xi,X2, then Xi + X2 
must correspond to Xi -\- X2 and X1X2 must correspond to XiX2- The notion of a 
*-algebra implements exactly the question posed above: every commutative *-algebra 
can be mapped to a set of random variables of the form £°°{!F) for some cr-algebra 
by using the spectral theorem and the following Lemma. The fact that sets of 
measurable functions can be characterized by algebras is well-known; in fact. Lemma 
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12.41 is just a "trivial" version of the monotone class theorem [721 Ch. I, Thm. 8]. 

Lemma 2.4. Let il = Then there is a one-to-one correspondence 

between (commutative) *-algehras of n x n diagonal matrices and £°° (J-)- spaces. 

Proof. We need two simple facts. First, let Xi, . . . , Xk be elements of a *-algebra 
of diagonal matrices. Then f{Xi, . . . , Xk) is also an element of the *-algebra for any 
function /. Clearly this is true if / is any (complex) polynomial. But this is sufficient, 
because given a finite number of points xi, . . . ,Xn G C*^ we can always find for any 
function / a polynomial / that coincides with / on xi, . . . , a;„. Second, we claim that 
for any there exists a finite set of disjoint sets Si, . . . , Sk & T such that IJj, Sk — ^ 
and such that any x G £°°{T) can be written as a; = J2i^iXSi ixs is the indicator 
function on S, Xi € C). To see this, note that there is only a finite possible number 
of disjoint partitions {Ti} of the latter being a finite set. {Si} is then the (unique) 
finest such partition that is a subset of J^. To show uniqueness: if {Si}, {S'^} C J- were 
two such partitions then {Si H S'j} d J-' would be a finer partition unless {Si} = {S^}. 

Let us now prove the Lemma. For a diagonal matrix X we define the map 
x{i) = L(X){i) = Xii, and similarly for a map a; : f2 — > C define the diagonal matrix 
Xii ~ L^^{x)ii ~ x{i). This gives a *-isomorphism l between the set of all diagonal 
matrices and the set £°°{^l) of all maps on fl. The claim of the Lemma is that t 
maps any *-subalgebra C ^ to £°°{J-) for some cr-algebra JF, and conversely that 
i~^(£°°(jF)) is a *-algebra for any !F. The latter is a straightforward consequence 
of the definitions, so let us prove the former statement. Fix the *-algebra ^ and 
define C = i(^). Let !F = ct{C}, and introduce the smallest il-partition {Si} C T as 
above. If xSi S C for every i, then we must have C — £°°{!F): C is an algebra and 
hence contains all random variables of the form X^i^iXSi- Let us thus assume that 
there is some i such that xSi ^ C\ as T = cr{C}, however, we must be able to find 
functions a;(i), . . . ,a;(„) G C such that x^^^{xi) n • • • nxj~^j(a;„) — Si. But if we choose 

f{yi, ■■■,yn) = X{yi=xi,...,y„=x„}{yi, ■ ■ ■,yn) then /(a;(i), . . . ,a;(„)) = xs, e C. Hence 
we have a contradiction, and the Lemma is proved. □ 

We are now ready to introduce the basic mathematical object studied in this 
article: a generalized or quantum probability space. Applying the spectral theorem to 
this structure gives a fundamental result which wc will use over and over. 

Definition 2.5 (Quantum probability space). The pair {£/,p), where £/ is 
a * -algebra of operators on a finite- dimensional Hilbert space H and p is a state on 
, is called a (finite- dimensional) quantum probability space. 

Theorem 2.6 (Spectral theorem for quantum probability spaces). Let 

, p) he a commutative quantum probability space. Then there is a probability space 
(n,.F,P) and a * -isomorphism l:'^ ^ £°°{T) such that p{X) = Ep(t(X)) VX e 

What is so general about a generalized or quantum probability space? It is the 
existence of many commutative subalgebras within (^,p). Theorem 12.61 does not 
apply directly to (^,p), as usually such a space will not be commutative. On the 
other hand, in a single realization we can choose to observe at most a commutative set 
of observables, which generate a commutative subalgebra ^ C £^ . The probability 
space ('^, p\cg) is commutative and is thus exactly equivalent to a classical probability 
space by Theorem 12.61 The noncommutative probability space (jz/, p) describes the 
statistics of all possible experiments — it is a collection of many incompatible classical 
probability models, each of which coincides with a commutative subalgebra of si/ . 
The experiment we choose to perform in one single realization determines which com- 
mutative subalgebra of , i.e. which classical probability space, is needed to describe 
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the random outcomes of that experiment. 

Remark 2.7. In section we determined the probabiHty measure P using a 
state vector ip ^ H; here we have replaced this notion by the state p : ^ ^ This 
is in fact a generahzation: the vector -0 corresponds to the state p{X) — {ip,Xip). 
In general we can always characterize a state p by a density matrix p as follows: 
p{X) = Tr[pX]. This follows directly from linearity of p. Positivity and normalization 
impose the additional conditions p > and Tr p — I. A state of the form p{X) = 
(ipjX^jj) = Tr[('0'i/'*)-'f] is known as a pure or vector state, whereas any other state is 
known as a mixed state. Both state vectors and density matrices are commonly used 
in the physics literature, while p : ^ C is usual in quantum probability. 

We will often speak of the *-algebra generated by a set of observables Xi, . . . , Xk- 
By this we mean the smallest *-algebra a\g{Xi, . . . , Xk} of operators on H that con- 
tains Xi, . . . , Xk- This notion plays exactly the same role as the cr-algebra generated 
by a set of random variables — indeed, it is straightforward to verify that if Xi , . . . , Xk 
commute, then t(alg{Xi, . . .,Xk}) = ^°°((7{i(Xi), . . . , t(Xfc)}). 

In the following we will need to construct quantum probability spaces from a 
number of independent smaller probability spaces. The analogous classical notion 
can be illustrated as follows: suppose (57^, J-^, P^), i = 1,2 are two independent copies 
of the probability space of a coin flip. Then (fii x ^2,^1 x ^2, Pi x P2) is the joint 
probability space for the two independent coin flips, on which any random variable 
f(uji,u!2) — /(wi) which depends only on the first coin and g{uJi,LU2) — 5(1-^2) which 
depends only on the second coin are independent (under the product measure Pi XP2). 

The analogous construction for quantum probability spaces uses the tensor prod- 
uct; that is, given quantum probability spaces {s^i^pi) defined on the Hilbert spaces 
Hi, we construct the joint space [s/i ® £^2iP\ ® P2) on the Hilbert space Hi (g) H2. 
Here Hi (g) H2 and s^i (g s^2 denote the usual vector and matrix tensor products, while 
Pi® P2 is defined by {pi® p2){A® B) — pi{A)p2{B) and is extended to all of jz/iiS)s^2 
by linearity. It is not difficult to verify that also this construction is consistent with 
the classical case, as £°°{!Fi x T2) = ® 1°°{T2)\ the *-isomorphism l obtained 

through the spectral theorem then maps observables of the form A® I to random 
variables of the form f{uji,uj2) — /(wi), and vice versa. 

For the time being this is all the general theory that we need. We now turn to 
the construction of the class of models that we will consider throughout the paper. 

2.3. Two-level systems. In this article all quantum systems will be built up 
(by means of tensor products) from two-level systems, i.e. quantum probability spaces 
with dim H = 2. Any observable of a two- level system is at most a two-state random 
variable, which makes the theory particularly simple. Nonetheless we will find a 
surprisingly rich structure in the resulting models. Moreover, it turns out that such 
models can approximate closely the behavior of realistic physical systems, see section 
1101 The classical counterpart of such models, i.e. the approximation of a continuous 
system by a sequence of coin flips, is well known e.g. in mathematical finance |261l75j. 

We already encountered a two-level system in section 12.11 As we will use these 
systems so often, however, it will be useful to fix some notation and to introduce some 
standard vectors and operators which we will use throughout. 

Let ^ denote the *-algebra of all 2 x 2 complex matrices acting on the Hilbert 
space H = C^. This algebra, together with suitable choices of states, will be the 
building block in our model of an atom in interaction with an electromagnetic field: 
the atom will be built on ^ (see Example I2.8p . while the field will be modeled by a 
tensor product of ./#'s (see section [^?5|) . 
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The canonical two-level system is described by the quantum probability space 
(^jP), where p is some state on Implicit in this description is the choice of a 
standard basis (hence the name canonical) which allows us to represent linear opera- 
tors by the more convenient matrices; there is no loss of generality in doing this, and 
it will make our life much easier. We introduce the standard definitions: 

o). --'^^(o J). 

$ is called the vacuum vector or ground state vector for reasons which will become 
clear shortly. Note that cr_<i? = 0, cr^ = 0, and ct_ct+ -|- ct+CT- = /. cr+ and cr_ are 
convenient operators, as any matrix M G ^ can be written as a linear combination 
of cr+<T_, (T_|-, (T_, /. We will also sometimes use the Pauli matrix 

a, = (T+cr_ - = . (2.2) 

Example 2.8. (Two- level atom). In this example we are going to interpret the 
quantum probability space (^, p) as the set of observables describing a single atom. 
Introduce the self-adjoint operator H E ^ by 

ii=—-^-\l (2-3) 

where ti is Planck's constant (a fundamental constant of nature), is a parameter 
representing a frequency (the so-called atomic frequency), and CTz is the Pauli matrix 
given by Eq. (|2.2p . As an observable, the Hamiltonian H is interpreted as the energy 
of the atom (indeed, hioo has units of energy). If an experiment is done to measure 
it, then from Eq. (|2.3|) it is clear that the outcome is a random variable taking either 
the value hujo/2 or — fia;o/2 (note that H is already diagonal in the given basis). 

The fact that the atomic energy takes discrete values is a fundamental result in 
quantum mechanics to which the theory owes its name. No actual atom has only 
two energy levels, but in some cases atomic physics is well approximated by such 
a description. Other physical observables, such as the intrinsic angular momentum 
(spin) of an electron, are exactly described by a two-level system. Throughout this 
article we will use the two-level system as a prototypical example of an atom. 

We have not yet discussed the state p; the state depends on the way the atom was 
prepared (e.g. is the atom at room temperature, has it been cooled to absolute zero, 
or perhaps it was excited by a laser?) A particularly interesting state, the ground 
state, is the one in which the energy takes its lowest possible value —huJo/2 with 
unit probability. This situation describes an atom that has been cooled to absolute 
zero temperature (a good approximation for laser-cooled atoms). It is not difficult 
to verify that the ground state is defined by p{X) = X^) where $ is the ground 
state vector. In general we can choose any state for the atom, i.e. we set 

p{X) = Tr[pX] with p = f^" , X G^, 

\P21 P22 J 

where p is a density matrix [p > 0, Tip ~ 1). Then measurement of H yields the 
outcome —hujQ/2 with probability p22 and hujQ/2 with probability pn. Note that the 
expectation value of the energy H is given by p{H) ~ Tr[/3_ff] ~ ^fkJo{pii — ^22)- 
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Fig. 2.1. Cartoon of a discretized model for the electromagnetic field. The field (e.g. a laser 
beam) propagates to the left until it hits a photodetector. When we discretize time with sufficiently 
small At = then to good approximation the photodetector will detect either zero or one photons in 
each time step, and moreover the number of photons detected in each time interval are independent 
random variables. The observable AA(i) represents the number of photons detected in the time 
interval i. A possible sample path for the stochastic process t(AA{j)) is displayed in red. 

Apart from the energy H the two-level atom may possess an electric dipole mo- 
ment, described by the vector observable /i. Its three components are given by 



where a^, OLy and az are complex parameters (depending on which specific atom we 
are considering). Note that a G {a;, y, z} does not commute with H . Therefore ^.a 
and H cannot be simultaneously diagonalized by the spectral theorem, that is they 
cannot be represented on the same probability space. It is impossible to devise an 
experiment that measures both observables simultaneously. Note further that ^x.y,z 
do not commute with each other unless ctx^y^z have the same phase, i.e. it is in general 
not even possible to measure all three components of ft in one realization. 

2.4. The discrete field. Having obtained a simple model for an atom, we pro- 
ceed by introducing a simple model for the electromagnetic field, e.g. light emitted 
from a laser. In the next section we will couple the atom and the field, which will 
allow us to describe measurements of light emitted from or scattered by the atom. 
Before we can do this, however, we must model the field in absence of the atom. 

The model we have in mind is illustrated in Figure 12.11 Imagine a laser beam 
which we are measuring with a photodetector. The beam propagates towards the 
photodetector at the speed of light, so at each time t we will measure a different part 
of the beam. We wish to describe the measurements made by the photodetector at 
each time t. To simplify matters we will discretize time into slices. If these slices are 
sufficiently small (as compared to the intensity of the laser), then to good approxima- 
tion the photodetector will measure no more than one photon in each slice. Moreover, 
the measurement in each slice will be independent from the other slices, as we are 
measuring independent parts of the field at each time step. Hence we will model 
the photodetection of the field by a collection of independent {0, l}-valued random 
variables, one for each time slice, corresponding to whether a photon was (1) or was 
not (0) detected in each slice. 

The model being a quantum mechanical model, the above description is not suf- 
ficient. The events detected by a photodetector are classical events, corresponding to 
an entire sample path of observations measured in a single realization, and hence are 
necessarily described by a commutative algebra. Other observations of the field are 
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possible, using a different detection setup, which are incompatible with direct pho- 
todetection (we will encounter an example below). Hence we need a noncommutative 
model. Still, the considerations above are very suggestive: we will build our model 
from a collection of independent two-level systems, one for each time slice. 

Fix an integer k and finite time interval [0, T], T > 0. Define A = yjT/k, i.e. the 
interval [0, T] is divided into k parts of equal length A^. To every time slice 1, . . . , fc we 
want to associate an independent copy of the canonical two-level system p). But 
we have already seen how to do this: we must use tensor products, i.e. we introduce 
the quantum probability space 

k times k times 

For any X e and integer 1 < « < fc, we define an element Xi ^ Wthy 

X, ®X®I'^^^~'\ 

Xi is the observable X in the time slice i. Note the fundamental property [Xi, Yj] = 
y X,Y G i ^ j- measurements on different slices of the field are always compatible. 
For future reference, we also define the algebra of observables in time slice i: 



Ji, = alg{X, : X e 
Discrete noises. For 1 < z < fc we define the standard increments 



(2-4) 



AA*(i) = A(fT+)„ A<(z) = A^/, 

and consequently we introduce the discrete noise processes A, A* , k and the discrete 
time process t (with the convention A(0) = ^*(0) = A(0) = t(0) = 0) as follows: 



A(Z) = ^A^(z), A(/)^^AA(0 

I I 

A*(0 = 5]A^*(z), t(0 = ^At(z) 



i=l 1=1 

The particular scaling of these expressions with A has been selected so that the def- 
initions make sense in the continuous limit A ^ 0; see section [TOl For now, let us 
furnish these objects with a suitable interpretation. t{l) = IX^I is the easiest: this 
observable takes the value IX^ with probability one. We interpret t{l) as the time that 
has elapsed after time slice I. Next, let us investigate A(/). 

Photon counting. For the observable a^a^ E it is easily verified that the 
spectral theorem maps to a random variable L{a^a-) that takes the value zero or one; 
indeed, the matrix corresponding to ct+CT- (calculated using Eq. (|2.ip ) is given by 



1 




which is already in diagonal form. AA(i), « = 1, . . . . fc is a set of independent ob- 
servables taking the value zero or one. We interpret AA(i) as the number of photons 
observed by a photodetector in the time slice i (as illustrated in Fig. 12. ip . 
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The probability of observing a photon in the time shce i depends on the state 
p; in the presence of a hght source (a laser or a light bulb) this probability will be 
nonzero. If we turn off the light source, or if there was no light source in the first 
place, then the probability of measuring a photon will be zero in every time slice. 
As we discussed in the context of a two-level atom, this situation is described by the 
state 4>{X) = (i.e. pk = <p®'^) — it is for this reason that this state is called 

the vacuum state. Even in the absence of an external light source an atom can still 
interact with the vacuum: e.g. an excited atom can spontaneously emit a photon into 
the vacuum, an event which we will want to measure using our photodetector. For 
concreteness, we will always work in the following with the vacuum state (f)®*' . 

Now recall that AA(?;), i — l,...,k all commute with each other, and hence 
generate a commutative *-algebra "^k- The spectral theorem maps % to a classical 
space {flk,J'k), which is simply the probability space corresponding to k independent 
coin flips (a "binomial model" in the terminology of 75]). The classical stochas- 
tic process yi = l{AA{1)) is then precisely the signal we observe coming from the 
photodetector. By applying the spectral theorem to the commutative subalgebras 
'rfi = alg{AA(i) : i — 1, . . . ,1}, '^i C %+i C • • • C ^fe, we obtain an increasing family 
of u-algebras J-) C J^i+i <Z ■ ■ ■ <Z Tk, the filtration generated by the stochastic process 
yi. As in the classical theory, the quantum filtration {"^i} allows us to keep track of 
what information is available based on observations up to and including time slice I. 

It remains to consider A(Z) € ^i: it is evident from the definition that this ob- 
servable represents the number of photons observed up to and including time slice I. 
Hence the number process {A(Z)};=i....^fc, which is clearly commutative as A(Z) g % 
V^, maps to a classical counting process t(A(Z)). All of this is not so interesting when 
we are working under the measure P/c induced by the vacuum (j)®^ on {^kT^^k)'- under 
this measure any sample path of /,(A(Z)) that is not identically zero has zero proba- 
bility. Once we have introduced the interaction with an atom, however, the number 
process will allow us to count the number of photons emitted by the atom. 

Homodyne detection. The processes A{1) and A*{1) are not self-adjoint, but 
X^p{l) — e'-'^A{l) + e^*'^A*(Z) is an observable for any phase e [0, tt). Moreover, 
{X(p(Z)};=i^....fc is a commutative process for fixed Lp: to see this, note that /S.Xif,{l) ~ 
Xip{l) — X^(l — 1) = A(e*'^(T_ -|- e~''''(T+)i commute with each other for different I, so 
the same must hold for X^p{l). Hence, as above, we can use the spectral theorem to 
map 'lo^ = &\g{X^{l) : / = 1, . . . , fc} to a classical space (fi^, J^^), and similarly we 
obtain the filtration {J'f}i=i^,,,^k- Beware, however: {Xi^} does not commute with 
{A}, nor does it commute with {X^p'} for Lp' ^ tp. To observe each of these processes 
we need fundamentally different detectors, so that in any given realization at most 
one of these processes is available to us. 

An optical detector that measures one of the processes X^{1) is known as a ho- 
modyne detection setup, a term borrowed from radio engineering. What is measured 
in such a setup is usually interpreted by considering the electromagnetic field to be a 
wave at a certain frequency, as in classical electromagnetism. The homodyne detector 
then measures the component of this wave which is in phase with a certain reference 
wave, known as the local oscillator, and the parameter ip determines the phase of the 
local oscillator with respect to an external reference. The interpretation of the light 
in terms of waves (homodyne) or particles (photodetection) are not at odds, because 
the corresponding measurements {A^,^} and {A} do not commute. We will not dwell 
further on the physical interpretation of these detection schemes; suffice it to say that 
both homodyne detection and photodetection are extremely common techniques, ei- 
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ther of which may be more convenient in any particular situation. The interested 
reader is referred to the quantum optics hterature [80] for further information (see 
also for a more mathematical perspective). 

For concreteness we will always consider the case ip = whenever we discuss 
homodyne detection in the following, i.e. we will consider the observation process 
X{1) = A{1) + A*{1). Let us investigate this process a little more closely. First, 
consider the two-level system observable cr+ + (T_ : 




P±i are projections, so we can read off that (t+ + cr_ has eigenvalues ±1. Moreover, 
note that ($, P±i^) = 1/2. Hence evidently the spectral theorem maps ct-i- + (T- to a 
random variable taking the values ±1 each with probability 1/2 in the vacuum. It fol- 
lows immediately that l{AX{1)), Z = 1, . . . , fc are independent random variables taking 
values ±A with equal probability. Unlike in the case of photodetection, the homodyne 
detector gives a noisy signal even in the vacuum: a manifestation of vacuum fluctua- 
tions, as this is called in quantum optics. The process l{X{1)) = X]i=i ^{^X^ (."i-)) j the 
integrated photocurrent, is evidently a symmetric random walk on the grid AZ. It is 
not difficult to see that i{X {[t / \)) converges in law to a Wiener process as A ^ 
(recall t{l) = iX^, so l{t) ~ [t/A^J)— see sectionHUl 

Remark 2.9. The reader might wonder at this point why we have chosen the 
notation i(AA(l)) or i,{AX{l)) for the observations — why not dispose of the A's? The 
reason is that neither AX{lt/X'^\) nor AA([t/A2j) have well-defined hmits as A ^ 0; 
for example, /,(AX([t/A^J )) would have to converge to white noise, which is not a 
mathematically well-defined object (at least in the sense of stochastic processes). The 
Wiener process, on the other hand, has a rigorous definition. We use the convention 
that processes such as X([</A^J) have well-defined limits, whereas AX{[t/X^\) will 
"converge" to dXf which is only meaningful under the integral sign. 

Remark 2.10. It should be emphasized that the model we have introduced 
here is rather crude; certainly it was not derived from physical principles! Its main 
physical justification is that it converges in the limit A — > to a model which is 
obtained from physical principles, so that for small A we can consider it to be a good 
approximation — see the references in section 1101 The same disclaimer holds for what 
we will discuss in the remainder of section [2] 

2.5. Repeated interaction. Now that we have descriptions of an atom and of 
the electromagnetic field, it is time to describe how these objects interact with each 
other. The atom can emit photons into the vacuum, which we subsequently detect 
using a photodetector; alternatively, we can measure using a homodyne detection 
setup how the presence of the atom perturbs the usual vacuum fiuctuations. The 
remainder of the article is centered around the theme: how can we put such measure- 
ments to good use? The development of filtering theory allows us to estimate atomic 
observables based on the measurements in the field, and we can subsequently use this 
information to control the atom by means of feedback control. 

The way the interaction works is illustrated in Figure 12.21 As before, the field 
propagates towards the detector, which measures sequentially time slices of the field. 
Now, however, we place the atom right in front of the detector, so that the current 
time slice interacts with the atom prior to being detected. This sort of model is called 
a repeated interaction: the atom is fixed, and interacts identically with each time slice 
of the field exactly once before it is detected. Before we can describe the repeated 
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i-1 2 3 4 i=2 ~ 3 4 5 

Fig. 2.2. Cartoon of the discretized atom-field interaction. In each time step i, the atom and 
the field at time slice i interact through a unitary transformation Mi. The time slice of the field 
is subsequently measured (using a photodetector or homodyne detection), which gives rise to the 
observation AY{i). The field then propagates to the left by one time slice, and the procedure repeats. 



interaction in detail, we first have to describe how a single interaction is modeled 
in quantum mechanics. Though this article is mostly about quantum probability, we 
need a little bit of quantum dynamics in order to model physical systems. 

Interactions in quantum mechanics. Suppose we have two independent two- 
level atoms, i.e. we are working on the quantum probability space p ® p). 
Then, for example, the observable p\ = p^ ® I corresponds to the electric dipole 
moment of the first atom, and is independent from any observable I®X of the second 
atom. If the two atoms subsequently interact, however, then the dipole moment of 
the first atom would likely change; i.e. the dipole moment after the interaction is 
described by a different observable than before the interaction, and will likely 
be correlated with some observable of the second atom (/ ® XY^^'^'^ . 

In quantum mechanics, every interaction is described by a unitary operator. In the 
previous example, the interaction between the two atoms corresponds to a suitably 
chosen unitary map U : ® €? ® . This means that any observable 

X G ^ (g) ^ transforms to U*XU S ^ (g) ^ after the two atoms have interacted. 
In particular, /iij.''^^''^'' = U* p^XJ will generally be different than p^, and is likely to be 
correlated with U* {I ® X)U for some X. Note that unitary rotation does not change 
the spectrum of an observable; i.e., l{pI.) and take the same values, but 

they are defined on different probability spaces with different probability measures. 

Recall that any unitary matrix U can be written &sU — e^*^ for some self-adjoint 
matrix B. Hence we can express any unitary transformation of ^ ® ^ by 

U — exp (— « {Li ® a+a^ + L2 ® (t+ + L*2 ® + ® I}) , 

where ii.2,3 G ^ and Li,Lj, are self-adjoint. 

Remark 2.11. We have described interactions by changing the observables 
corresponding to a particular physical quantity. This is how dynamics is usually 
represented in probability theory. For example, the classical discrete-time Markov 
chain X„ = /(X„_i,^„) where ^„ is i.i.d. noise is precisely of this form: the random 
variables X„ and X„_i represent the same physical quantity X, but the interaction 
with the noise ^ means we have to represent X by different random variables Xn 
at different times. In quantum mechanics, this is called the Heisenberg picture. In 
the Schrddinger picture, any physical quantity is always represented by the same 
observable, but the underlying state is changed by the interaction. The two pictures 
are in some sense dual to each other, as is discussed in any quantum mechanics 
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textbook (e.g. j68j). However, the Schrodinger picture is somewhat unnatural if one 
wants to work with (quantum) stochastic processes. In this article we will always 
work in the Heisenberg picture. 

Atom-field interaction. We now consider the case that a two-level atom is 
coupled to the vacuum electromagnetic field, that is, we work on the quantum prob- 
ability space {J( (g) P) with P = p (g) . The subspace {J(, p) C (^ (g) ^fe, P) is 
known as the initial system (p is called the initial state), because observables of the 
form X g) / in ^ (g Wk represent the atom's physical quantities at the initial time 
I = 0, i.e., before the atom has interacted with the field. The interaction requires us 
to modify these observables at each time Z > 0, as we will now describe. 

Let us begin by considering what happens to the system after one time step. At 
this time (I = 1) the atom has interacted only with the first time slice of the field, i.e. 
any atomic observable X®I aX time / = will evolve to ji{X) = M(l)*(Ag)/)M(l), 
where M{1) is the unitary interaction of the form 

M(l) = e-«{-f'iAA(l)+L2AA*(l)+L;AA(l)+L3®At(l)} 

In the next time step, we want the atom to interact in an identical fashion with the 
second time slice, so we will define i2{X) = M{2)*ji{X)M{2) with 

M(2) = g-«{ii(ii)AA(2)+ii(L2)AA-(2)+ji(L;)AA(2)+ii(L3)At(2)}^ 

Note that we do not only need to change the time slice that we are interacting with, 
but we also need to propagate the atomic observables that define the interaction. 

Remark 2.12. Once again it is useful to have in mind the example of the 
classical Markov chain X„ = /(A„_i,^„). Suppose X„ takes values in M'"'; if we 
are looking e.g. at the first component X^, then X^ — f{X^_i,X^_i, . . . ,X!^_i,£_n). 
Hence the transformation f{-,X!^_i, . . . ,X^_i,^n) of X^_i depends not only on the 
noise "time slice" but also on other "observables" of the system A^_;^, . . . ,X!^_^ 
which must be evaluated at time n—1. Similarly our quantum mechanical interaction 
M depends at each time step I on the Ith time slice in the field as well as on the 
atomic observables in the previous time step I — 1. 

Proceeding in a similar fashion, we find that any atomic observable X e ^ has 
evolved to the observable ji{X) = U{l)*{X g) I)U{1) £ Ji( ®Wk sX time I, where 

U{1) = n^=iM(z) = M(1)M(2) • • • M{1), f/(0) = /, (2.5) 
is a repeated interaction with the one time step interaction 

M{1) = (,-i{]i-i(Li)AHl)+Ji-i(L2)^A'(l)+]i_i{Ll)^A{l)+ji_i{L3)At{l)} _ ^2.6) 

The map defines the time evolution or flow of the atomic observables. The 
connection with dynamical systems theory can be made explicit: if we define U{i,l) = 
Mii + l)M{i + 2)...M{l), i < I, and J^,i{X) = U{iJ)*XU{iJ) for A d Ji ®Wk, 
then Ji,;(-) is a two-parameter group of transformations of the algebra ^ (g #fc (i.e. 
for i < ; < r Ji,riJi,i{-)) = Ji,r{-), JiA'Kli )) = JlA^)^ etc.); thus we have truly 
defined a discrete time dynamical system in the dynamical systems sense. We will 
not need this in the following, however. 

There is a different representation of the repeated interaction unitary U (/) which 
is often more convenient. Introduce the unitaries 



j^r ^ g-i{LiAA(0+L2AA*(0+L;AA(0+L3At(0} 



(2.7) 



A DISCRETE INVITATION TO QUANTUM FILTERING AND FEEDBACK CONTROL 17 



which, in contrast to M{1), depend on the initial observables ^1,2,3 rather than the 
time-evolved observable ji_i(Li^2, 3)- From j;_i(ii.2,3) = U{l~l)*{Li^2.3'^I)U{l~l), 
it follows immediately that M{1) = U{1 - l)*MiU{l ~ 1). But then 

Uil) = U{1 - l)M{l) = U{1 - l)U{l - l)*MiU{l - 1) = MiU{l - 1), 

so we obtain the expression 

U{1) = n^=iM, = A//Mz_i • • • M2M1. 

Hence we can use the A/;'s rather than M(l)'s if we reverse the time ordering in the 
definition of U{1). This is often convenient because Mi is only a function of the initial 
system and of the ^th time slice, whereas M{1) depends on the entire history of the 
field up to time slice I. 

The choice of L12.3 € ^ determines the nature of the physical interaction be- 
tween the atom and the field. For different atoms or for different experimental scenar- 
ios this interaction can take very different forms; nonetheless the theory that we are 
about to develop can be set up in quite a general way for a large class of interactions. 
In section 12.61 we will introduce two specific examples of physical interactions which 
we will use throughout to illustrate the general theory. 

Observations. Up to this point we have concentrated on the time evolution 
ji{X) of an atomic observable X ^ I; however, the observables corresponding to time 
slices of the field are also modified by the interaction with the atom in a completely 
analogous way. Recall, for example, that before the interaction the number of photons 
in time slice / is represented by the observable J® AA(/). After the interaction this 
quantity is represented by AY^{1) ^ U{l)*{I <S> AA{l))U{l), so that l{AY^{1)) is 
precisely what we observe at time I if a photodctcctor is used to measure the field. 
Similarly, homodyne detection measures AY^{1) = U (l)* (I AA{1) + 1 <^ AA* {l))U (1) 
at time /. We wiU use the generic notation AY{1) = U (l)* AZ {l)U (l) , where AY = 
AY^ ifAZ = I(g)AA and AY = AY^ if AZ = / ® {AA + AA*). 

The first question we should ask ourselves is: does this procedure make sense? 
After all, in the laboratory we measure classical stochastic processes, so that in order 
for our theory to be consistent the random variables l{AY{1)), I = 1, . . . ,k must define 
a classical stochastic process on some fixed probability space (fifc, .Ffc, P/c) — in other 
words, = alg{Ay(/) : I ~ 1, . . . , k} must be a commutative algebra. Let us verify 
that this is indeed the case. The basic insight we need is the following observation: 
[Ml, AZ{i)] = ioT I > i (this is easily verified by inspection), so that for j > I 

= (fi-i^^*) ^^(0 (J\!.=lM^^ = U{irAZ{l)U{l). 

This is not unexpected: it is simply a statement of the fact that the Ith time slice 
of the field interacts exactly once with the atom, viz. at the Ith time step (note that 
U{j)*AZ{l)U{j) = AZ[l) for j < I), for any /. But then for j > I 

[AY{l),AY{j)] = [U{irAZil)U{l),U{jrAZ{j)U{j)] 

= [UijrAZil)Uij),UijrAZij)Uij)] = UijnAZil),AZij)]Uij) = 0, 



18 



BOUTEN, VAN HANDEL, AND JAMES 



SO that clearly ^ is commutative. This self-nondemolition property guarantees that 
the observations can indeed be interpreted as a classical stochastic process through the 
spectral theorem, a crucial requirement for the physical interpretation of the theory. 

Previously we introduced the filtration '^i = alg{AZ(i) : i = 1, . . . ,1}, I = 1, . . . ,k 
corresponding to the photodetection or homodyne detection of the field in absence of 
interaction with the atom. Similarly we now define the filtration iVi = alg{Ay(i) : i = 
1, . . . , I — 1, . . . , fc which represents the information contained in observations of 
the field (after it has interacted with the atom) up to the Ith time step. Applying the 
spectral theorem to Wk, we obtain a classical probability space with the corresponding 
observation process Ay; — l{AY{1)) and filtration — i°°{yi) (note that we then 
have 3^; = (T{Ayi, . . . , Ayi}). The following fact will be useful later on: 

m = uurmij) - {uurxuij) -.xe^} vj > i. 

The proof is identical to the proof of the self-nondemolition property. 

We conclude this section with the demonstration of an important property of the 
models under investigation, the nondemolition property, to which we have already 
alluded in section [TT] Choose any X G ^ . Then for any I < i 

[AY{l),j,{X)] = [UiirAZil)U{l),U{^nx®I)U{i)] 

= [U{iyAZ{l)U{i),U{iyiX ® I)U{i)] ^U{iy[AZ{l),X ® I]U{i) =0. 

Evidently any atomic observable in time step i commutes with the entire history of 
observations up to that time. Thus in principle we could decide at any time to stop 
observing the field and we would still be able to measure any atomic observable: the 
joint distribution of AY{1), . . . , AY{1) and ji{X) is defined for any self-adjoint X e 
as alg{AF(l), . . . , AY{1), ji{X)} is a commutative *-algebra, despite the fact that 
ji{Xi) and ji{X2) need not commute for Xi X2- This enables us to meaningfully 
define the conditional expectation of ji{X) with respect to the observations AY(i), 
i = 1, . . . J for any X € which we will do in section[3l The nondemolition property 
provides us with a sensible way of estimating a whole set observables, despite the fact 
that these do not commute with each other, as every observable in this set separately 
commutes with the observation history on which we are basing the estimate. If the 
latter were not the case the estimates would have no physical relevance: one would not 
be able to devise an experiment, even in principle, which could verify the predictions 
of such an estimator. Similarly, it would be meaningless to try to control atomic 
observables which do not commute with the observations on which we have based our 
feedback control law. The nondemolition property avoids such problems and is thus 
of fundamental importance for quantum filtering and feedback control. 

2.6. Examples. The interaction matrices ^1,2,3 determine the nature of the 
physical interaction between the atom and the field. Though we will set up the theory 
in the following sections for arbitrary ^1^2, 3, we will repeatedly demonstrate the theory 
using the following two examples. The examples illustrate two common experimental 
scenarios — spontaneous emission and dispersive interaction — and we will be able to 
study these examples explicitly through numerical simulations. 

Spontaneous emission. The spontaneous emission scenario is obtained when 
an atom sits undisturbed in the vacuum at zero temperature. If the energy of the 
atom is minimal ¥{H) = —hLUQ/2, then it will remain this way for all time. Otherwise, 
the energy of the atom decays to its minimal value. A photodetector measuring the 
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field would see exactly one emitted photon at some randomly distributed time: hence 
the name spontaneous emission. Wc will reproduce this behavior through numerical 
simulations once we have introduced the corresponding filtering equations. 

The spontaneous emission model is characterized by the interaction matrices 



where n is the spontaneous emission rate. For simplicity, wc will always set 2k = 1. 
This is a discretized version of the well known Wigner-Weisskopf model. 

Remark 2.13. In principle, we could calculate the possible observation sample 
paths and their probabilities already at this point. After all, we would only need 
to simultaneously diagonalize the matrices AF(i), i = l,...,fc obtained using this 
explicit choice for ^1,2,3- We will have a much more efficient way of calculating such 
sample paths, however, once we have introduced the filtering equations: these will 
allow us to simulate typical sample paths using a Monte Carlo method. Hence we 
postpone numerical investigation of our examples until then. 

Dispersive interaction. In the spontaneous emission scenario, the atom ra- 
diates its excess energy by emitting a photon with energy fkoo (i.e. with frequency 
loq) — this is precisely the difference between the two values that the energy H can 
take. It is possible to suppress the spontaneous emission significantly by "discour- 
aging" the atom to emit photons with frequency wq; for example, wc can place the 
atom in an optical cavity whose resonant frequency is far detuned from ujq. In this 
dispersive regime, the (effective) interaction between the atom and the field is rather 
different in nature. An atom which has energy +fta;o/2 will shift the phase of the 
output light by a small amount in one direction, while an atom with energy —fujJo/2 
will shift the phase by the same amount in the other direction. Such phase shifts can 
be observed using a homodyne detection setup, and the resulting photocurrent thus 
carries information on the energy of the atom. 

The dispersive interaction model is characterized by the interaction matrices 



where g is the interaction strength. For simplicity, we will always set (7 = 1. 

3. Conditional expectations and the filtering problem. Within the con- 
text of the repeated interaction models introduced in the previous section, we have 
established that the observation history up to time step I is compatible with any 
atomic observable at time step I. In this section we will show how to estimate atomic 
observables based on the observation history. Because only commuting observables are 
involved, the corresponding theory can simply be "lifted" from classical probability 
theory using the spectral theorem. This is precisely what wc will do. 

3.1. A brief reminder. We begin by briefly recalling how classical conditioning 
works. Let (fi, J^, P) be a probability space and let P, Q e be a pair of events. The 
conditional probability of P given Q is given by 



Li=L3 = 



L2 = iV^Ka-, 



Li = L3 = 0, 



L2 = i\fg<J. 



P(png) 
P(Q) 



provided P(Q) 7^ 0. One could interpret this quantity roughly as follows: generate a 
large number of samples distributed according to the measure P, but discard those 
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samples for which the event Q is false. Then P(P|Q) is the fraction of the remaining 
samples for which the event P is true. 

Now suppose f : Q. ^ {/i, . . . , /„} and g : Q. ^ {gi, . . . , g^} are measurable 
random variables that take a finite number of values. The conditional expectation of 
/ given g is the random variable defined by 

n m 

Ep(/|5)H=E^9-M9,)HE/'P(/"'(/*)l5"'(5,)), ^ef^- (3.1) 

j=l i=l 

In fact, we can consider Ep(/|5)(a;) to be a function of g: 

m 

Ep(/|5)M = F{g{cuj), F : g, ^ f,P{r\f,)\g-\g^)). 

1=1 

To interpret this quantity, again we generate a large number of samples but now we 
divide these up into disjoint subsets corresponding to the value taken by g. Now 
average / over each of these subsets. Then F{gj) is the average of / over the subset 
of samples on which g takes the value gj, so that Ep{f\g) = F{g) is simply the 
expectation of / given that we know g. Note that we did not define the quantity 
P(i-'IQ) for the case P((5) = 0; hence the expressions above for P(P|(5) are properly 
defined everywhere on O except on = Uj P(g-i(3j))=o ^ -^^ allow 

Ep (/!(;) to take arbitrary values on ^l^ and call any such random variable a version 
of the conditional expectation of / with respect to g. There is no loss in doing this: 
all versions of the conditional expectation coincide with unit probability. 

The use of these expressions is rather limited, as they do not extend beyond ran- 
dom variables g that take a finite number of values. One of Kolmogorov's fundamental 
insights in the axiomatization of probability theory was his abstract definition of the 
conditional expectation, which works in the most general settings. In this article we 
only use finite state random variables, but even so we find it significantly easier to 
use the abstract definition than these clumsy explicit expressions. Let us thus briefly 
recall the general setting, as can be found in any probability textbook (e.g. [82]). 

Let Q <Z J- hea. cr-algebra on which we want to condition — i.e., these are the events 
which we know to be true or false due to our observations. Let / be a J^-measurable 
random variable. Then any tj-measurable random variable / that satisfies 

Ep(/5) = Ep(/g) V ^-measurable (bounded) g (3.2) 

is a version of the conditional expectation / — Ep(/|tJ). The conditional expectation 
is by construction a function of what we have observed (it is ^-measurable), as it 
should be. It is easily verified that the explicit expression Eq. I|3.ip satisfies the 
abstract definition with Q = (7{g}. In fact, in the case where takes a finite number 
of values we can always express Ep(/|fJ) as in Eq. (|3.ip by constructing g specifically 
so that it generates Q: expressing g in terms of the smallest measurable partition of 
as in the proof of Lemma 12.41 we only need to make sure that g takes a different 
value on each set in that partition. Hence there always exists a / that satisfies the 
abstract definition, and it is not difficult to prove that all versions coincide with unit 
probability — the proof of the latter is the same as the one given below in the quantum 
case. Elementary properties of the conditional expectation are listed in Table [XT] 

We conclude by recalling the geometric interpretation of the conditional expec- 
tation. For simplicity, take 17 = {!,..., A^} and fix the cr-algebra Q <Z T. Clearly 
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Table 3.1 

Elementary properties of the conditional expectation in classical (and quantum) probability; see 
e.g. [82]. Equalities are always in the sense of versions, e.g. Ep (/i + Ep{/2 = Ep(/i + /2|5) 
means that any version o/Ep(/i|C7) + Ep(/2|C7) is a version o/Ep(/i + f2\Q). 



Property 


Description 


Linearity 
Positivity 

Invariance of expectation 
Module property 

Tower property 
Independence 


aEp(/i|g) +/3Ep(/2|g) = Ep(a/i +/3/2I6;). 
If / > then Ep(/|6) > a.s. 
Ep(Ep(/|e)) = Ep(/). 

If g is C7-measurablc, Ep(/g CJ) = gEp{f\Q) a.s. 
In particular, 'Ep{g\G) = (;Ep(l|5) = g a.s. 

liHcg, Ep(Ep(/|g)|w) = Ep(/|w). 

lfa{f,n} is independent of g, Ep (/|o-{W, 5}) = Ep(/|W) a.s. 



£^{J-') is a finite-dimensional linear space and £°°{G) C £°°{J-) is a linear subspace. 
Moreover {f,g)p = Ep{f*g) is a pre-inner product: it satisfies all the conditions 
of the inner product, except (/, /)p = ifi^ / = P-a.s. rather than iff /(w) = 
Voj G ri. Hence under (•, •)p, is a pre-inner product space. With this interpre- 

tation, it is evident from Eq. (|3.2p that Ep(/|C/) is simply the orthogonal projection 
of / G £°°{!F) onto i°°{G)- It follows from the projection theorem that Ep(/|5) is 
the a.s. unique random variable / e i°°{G) that minimizes the mean-square error 
{f — f, f — f)p = Ep(|/ — /P). This gives the conditional expectation, in addition to 
the probabilistic interpretation above, a firm statistical interpretation as the estimator 
that minimizes the mean square error. 

3.2. Quantum conditional expectation. Let {3§,P) be a commutative quan- 
tum probability space; by the spectral theorem, it is equivalent to some classical 
probability space (fi,^, P). Moreover, the *-subalgebra ^ C ^ defines through the 
spectral theorem a cr-subalgebra G C {as = 1°°{G) C £°°{J^)). Because the 
classical and quantum probability models are completely equivalent, we can simply 
lift the definition of the conditional expectation to the algebraic level: 

P(X|'^) = i-i(Ep(t(X)|g)) yX(^SS. 

This is nothing but a re-expression of the conditional expectation Ep(/|5) in the 
language of *-algebras. In fact, we can go one step further and directly translate 
the abstract definition Eq. (|3.2p into algebraic language: any X € ^ that satisfies 
f{XC) = P(XC) VC e is a version of ¥{XY^). We emphasize again: by the 
spectral theorem, observables X ^ SS arc just random variables and P(X|^) is just 
the ordinary conditional expectation. The spectral theorem is a powerful tool indeed! 

Usually we do not start with a commutative quantum probability space. Rather, 
we begin with a noncomniutative space (^, P) and choose a commutative subalgebra 
^ C ^ corresponding to the observations — keep in mind the example of the space 
p®(\)®^^ with the observations C ^®Wk at time step I. As we have seen, 
there could be many elements X ^ that commute with every element in but 
that do not commute with each other. The set ^' = {A e ^ : [A, C] = VC e is 
called the commutant of in si/ . It is easily verified that '^S" is always a *-algebra, but 
generally is not commutative. Nonetheless we can naturally extend our previous 
definition of P(-|'^^) from the commutative algebra ^ to the case ^ = 

Definition 3.1 (Quantum conditional expectation). Let (j2/,P) be a quan- 
tum probability space, and let '^S' be a commutative subalgebra of . Then the map 
P(-|'^) : ^ ^<a is called a version of the conditional expectation from onto ^ if 
P(P(yl|'r)C) = V{AC) yA e "T', C e "if. 
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How should we interpret this definition? Let A — A* ^ he an observable, 
and suppose we are interested in P(A|'^). The definition of P(j4|^) only involves the 
operators A and C € so for the purposes of evaluating P(A|^) we don't need 
the entire algebra we might as well confine ourselves to '^a = alg{^, C . C G 
■^l C "t^' — but 'i^A is commutative! Hence for any A — A* G 'rf', Definition 13.11 simplv 
reduces to the commutative definition we gave before. For A A*, we can always 
write A = P + iQ where P = {A + A*)/2 and Q = -i{A - A*)/2 are self-adjoint, 
and we define P(A|^) = F{P\'^) + iP(Q|^). We do this mostly for computational 
convenience (e.g., to ensure that the conditional expectation is linear); at the end of 
the day, only conditional expectations of observables have a meaningful interpretation. 

By now it should not come as a surprise that the conditional expectation can only 
be defined on the commutant . In section [3TT] we gave an "experimental procedure" 
for calculating the conditional expectation of a random variable /: generate a large 
ntmiber of samples from the joint distribution of / and the observations which we 
are conditioning on, separate the samples into subsets corresponding to the different 
observations, and average / within each of these subsets. But quantum mechanics 
tells us that it is fundamentally impossible to devise such an experiment when the 
observable A which we want to condition is incompatible with the observations! The 
commutant is thus precisely the right notion for a meaningful definition of the 
quantum conditional expectatior0. 

Having introduced an abstract definition, let's do a little ground work. 

Lemma 3.2. P(^|'^), A £ exists and is unique with probability one, i.e. any 
two versions P and Q of¥{A\'S') satisfy P((P - Q)*(P - Q)) = 0. 

Proof. We first introduce some notation. {F,G)p = P{F*G) defines a pre- inner 
product on thus turning it into a pre-inner product space. The associated semi- 
norm is denoted by ||i^||p = {{F,F)py/^. 

Existence. For any self-adjoint A^A*e'if', define "^a = alg{A, C : C e'lf}. The 
spectral theorem then provides a probability space {^a, J' a, Pa) and a *-isomorphism 
M : ^ e'^iJ'A), and we write mC^) = i°°{gA). For every such A, define P(A|'^) = 
t^^(a) where a is some version of Pa{i-a{A)\Qa). For non-sclfadjoint A G 'rf' , define 
¥{A\'^) = P(F|"r) + if'{Q\'^) with P ^ {A + A*)/2 and Q = - A*)/2. It is 
easily verified that the map P(-|'^^) thus constructed satisfies Definition 13.11 

Uniqueness with probability one. Let P and Q be two versions of P(A|'^). From 
Definition O it follows that (C, P - Q)p = for all C e But P - Q £ "T, so 
{P-Q,P-Q)p = \\P-Q\\^ = 0. □ 

Let us prove that the conditional expectation is the least mean square estimate. 

Lemma 3.3. P(A|'^) is the least mean square estimate of A given , i.e. for any 
C e'f we have \\A - P(A|'^^)||p < \\A - C\\p. 

Proof For any C € ^ we have \\A - C||^ = \\A - ¥{A\'^) + ¥{A\'^) - C||?. Now 
note that by Definition [SH we have (C, A - P(A|^))p = F{C*A) - F{C*¥{A\'^)) = 
for all C G ^, ^ e i.e., A-ViAl"^) is orthogonal to ^. In particular, yl - P^yll-r) 
is orthogonal to P(A|'r) - C e and we obtain \\A - C\\^ = \\A - F{A\'^)\\^ + 
||P(A|'^) - C\\l. The result follows immediately. □ 

The elementary properties listed in Table [3Tl and their proofs carry over directly 
to the quantum case. For example, let us prove linearity. It suffices to show that 



That is not to say that this is the only definition used in the quantum probabihty literature; in 
fact, a more general definition, of which our definition is a special case, is very common [78] . Such 
a definition is very different in spirit, however, and is not motivated by the probabilistic/statistical 
point of view needed for filtering and feedback control. 
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Z = aP(A|<^) + (3P{B\'^) satisfies the definition of P{aA + PB]'^), i.e. P(ZC) = 
F({aA + (3B)C) for all C G But this is immediate from the linearity of P and 
Definition 13. II We encourage the reader to verify the remaining properties. 

In Section [5] we will need to relate conditional expectations with respect to dif- 
ferent states to each other. This is done by the following Bayes-type formula. 

Lemma 3.4 (Bayes formula). Let (j2/,P) he a quantum probability space, and let 
be a commutative subalgebra oj . Let V be an element in such that V*V > 
and ¥{V*V) ^1. We can define a new state Q on by Q{X) = F(V*XV), and 



Proof. Let K be an element of For any X £ we can write 



F(F{V*XV\'^)K) = F{V*XKV) = Q{XK) = Q{Q{X\'^)K) = 

F{V*VQ{X\'^)K) ^ F{F{V*VQ{X\'^)K\'^)) = F{F{V*V\'^)Q{X\'^)K), 

and the result follows directly. □ 

3.3. Statement of the filtering problem. Let us return to the two-level sys- 
tem coupled to the discretized electromagnetic field, as described in section [2T5l Recall 
that we observe the commutative process AY{1), with either Y = or Y^ depending 
on what detector we choose to use, and that this process generates the commutative 
filtration '3^i. Suppose we have been observing Y up to and including time I. What 
can we infer about the atom? 

The nondemolition condition allows us to make sense of this question. In the 
language of this section, nondemolition simply means that ji{X) £ ^V/ for any I = 
0, . . . , k, X E ^ . Hence the conditional expectations 

^i(X) = F{2i{X)\%), Q<l<k, XEJi 

are well defined, and we could e.g. interpret i^i(X^ as the least mean square estimate 
of an atomic observable X at time step I given the observation history. 

Though it is evident that the conditional expectations i^iiX) are well defined, it 
is quite another matter to calculate them explicitly. In principle we have a recipe 
which we could follow to perform this calculation: find a basis in which ji(X) and 
Ay(i), i = 1,...,? are simultaneously diagonalized, apply an expression similar to 
Eq. p.ip . and then re-express the result in the original basis. Such a procedure would 
be extremely unpleasant, to say the least; clearly we are in need a better way of doing 
this. Finding a useful explicit representation of the conditional expectations 7ri(-) is 
known as the filtering problem. 

Sections [5] and [6] are dedicated to the solution of the filtering problem, which we 
will do in two different ways (and obtain two different expressions!) Both methods 
give rise to recursive equations, that is, we will find filtering equations of the form 
TTi{-) = /(7r;_i(-), Ay(/)). This is very convenient computationally as we don't need 
to remember the history of observations: to calculate the current estimates t^i{-), all 
we need to know are the corresponding estimates in the previous time step 7r/_i(-) 
and our observation in the current time step AY {I). In other words, we can update 
our estimates in real time using only the current observation in each time step. 
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3.4. A brief look ahead. The filtering problem is very interesting in itself; its 
solution provides us with an efhcient way of estimating atomic observables from the 
observation data. We have an ulterior motive, however, for putting such emphasis 
on this problem — filtering is a fundamental part of stochastic control with partial 
observations, and hence plays a central role in quantum feedback control. This is 
the subject of sections [THU and we postpone a detailed discussion until then. Let us 
briefly sketch, however, why we might expect filtering to enter the picture. 

In order to design a controller, we first have to specify what goal we are trying to 
achieve. Control goals are generally expressed in terms of expectations over atomic 
observables. For example, we might want to maximize the expected a;-dipole moment 
'^{jkiPx)) the terminal time t{k) — T, or we might try to steer the dipole moment 
to a specific value a (e.g. by minimizing P(jfc((/ia; — a)^))). If we are not sure that we 
are going to finish running the control through time T, perhaps we would prefer to 
maximize the time-average expected dipole moment ^ X)f=i ^Uiil^x)) — etc. etc. The 
problem is that our control can depend on the entire observation history; the optimal 
control can have quite a complicated dependence on the observation history! 

The key observation we need is that by an elementary property of the conditional 
expectation, P{ji{X)) = F{Tri{X)). Rather than expressing our control goal directly 
in terms of the quantum model, we are free to express the same goal in terms of the 
filter only, as though it were the filter we wanted to control rather than the atom 
itself. It will turn out that optimal controls have a very simple dependence on the 
filter: the control in time step I can be expressed as some (deterministic) function of 
the filter at time / — 1. The filter 7r;(-) is called an information state for the control 
problem, as 7r;(-) contains all the information we need to compute the feedback law. 
This is discussed in detail in section [5?^ 

The fact that optimal controls separate into a filtering step and a simple feed- 
back step is known as the separation principle. The complicated dependence of the 
control on the observation history is completely absorbed by the filter (which is by 
construction only a function of the observations). By solving the filtering problem in a 
recursive manner, we obtain an efficient way of calculating the optimal controls. This 
amply motivates, beside its intrinsic interest, our interest in the filtering problem. 

4. Discrete quantum stochastic calculus. In the previous sections we have 
provided an introduction to quantum probability, introduced a toy model for the 
interaction of an atom with the electromagnetic field, and introduced conditional 
expectations and the filtering problem in the quantum setting. At this point we 
change gears, and concentrate on developing machinery that will help us to actually 
solve the filtering problem (in two different ways) and ultimately to treat control 
problems associated with our model. 

We begin in this section by rewriting the repeated interactions of section 12.51 as 
difference equations, and we develop a "stochastic calculus" to manipulate such equa- 
tions. In the discrete setting this calculus is not a very deep result: it is merely a 
matrix multiplication table, as can be seen from the proof below. In continuous time, 
however, difference equations become (quantum stochastic) differential equations and 
the stochastic calculus is an indispensable tool. While in continuous time linear alge- 
bra is replaced by the much more difficult functional analysis, the quantum stochastic 
calculus nonetheless enables one to deal with the corresponding objects using simple 
algebraic manipulations. Bearing in mind our goal to mirror as closely as possible the 
continuous theory, we will consistently use the discrete calculus in what follows. 
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4.1. The discrete quantum Ito calculus. Let us fix once and for all the 
quantum probability space (8) #fe,P) with P — p (g) cj)^'', which models an atom 
together with the electromagnetic field as in section 12.51 The algebra ^ (g) #fe has a 
natural noncommutative filtration {S8i\i=i....k defined by 

i^o = ® z®'' c ^1 = ^ ^ c • • • c = ^ ® :rfe. 

The noncommutative algebra contains all observables that are a function of the 
atom and the time slices of the field up to and including slice I. 

A quantum process is a map from {0, 1, . . . , A:} to ^®y^k- A quantum process L is 
called adapted if L{i) E S^i for every < i < fc. It is called predictable if L(i) € 
for every 1 < i < fc. For a quantum process L we define AL as AL(i) — L{i) — L{i — 1) 
for 1 < i < fc. These definitions are similar to those used in classical probability [52] ■ 

Definition 4.1. (Discrete quantum stochastic integral) Let L be a pre- 
dictable quantum process and let M be one of the processes A, A* , A or t. The 
transform of M by L, denoted L ■ M is the adapted quantum process defined by 

I 

L- M{l) = ^L{i)AM{i), l<l<k, L-M(0) = 0. 

1=1 

L ■ M is also called the discrete quantum stochastic integral of L with respect to M . 

Note that for any discrete quantum stochastic integral as above, we can equiva- 
lently write X = L-M -^=> AX = LAM. Note also that L{i) and AM{i) commute 
by construction, so that we can always exchange the order of the (predictable) in- 
tegrands and the increments AA, AA, AA* and At. We can now state the main 
theorem of our "stochastic calculus" . 

Theorem 4.2. (Discrete quantum Ito rule) Let Li and L2 be predictable 
quantum processes, Mi, M2 € {A, A*, A, t}, and let X and Y be the transforms Li-Mi 
and L2 ■ M2, respectively. Then 

A{XY){1) = X{1 - 1) Ar(/) + AX{1) Y{1 - 1) + AX{1) AY{1), 

where AX{1) AY{1) should be evaluated according to the discrete quantum Ltd table: 



AX\AY 


AA 


AA 


A^* 


At 


AA 





AA 


At - X'AA 


X'^AA 


AA 





AA 


AA* 


A^AA 


AA* 


A^AA 








A^A^l* 


At 


X^AA 


A^AA 


X^AA* 


X^At 



For example, if Mi = A and M2 ^ A*, then AX AY = LiLaAt - A^LiLzAA. 

Proof. The quantum Ito rule follows from writing out the definition of A{XY): 

A{XY){1) = X{1) Y{1) - X{1 - 1) Y{1 - 1) 

= {X{1 - 1) + AX{l)){Y{l - 1) + AY{1)) - X{1 - 1) Y{1 - 1) 
= X{1 - 1) Ar(/) + AA(/) Y{1 - 1) + AX{1) AY{1). 

The Ito table follows directly from the definition of the increments A A, A A, AA* , At 
in section [2.41 bv explicit matrix multiplication. □ 
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The following lemma shows that transforms of A, A* and A have zero expectation. 
It is a fancy way of saying that (t>{(J-) = = ^(o'+f-) — and (/>(/) = 1. 

Lemma 4.3. For all predictable processes L we have F{L{l)At{l)) — P(i(/)) 
and V{L{l)AA{l)) = F{L{l)AA*{l)) = P(L(/)AA(/)) = for any 1 < I < k. 

Proof. Note that (/)((T-) = 't'{<^+) = '/'(o'+o'-) — and — 1. Since L is 

predictable, L{i) and AA{i) are independent under P and we have for < i < fc 

f{L{i)AA{i)) = AP(i(i)) 0(cr_) = 0. 

Similar reasoning holds for A* , A and t. □ 

4.2. Quantum stochastic difference equations. We are now going to rewrite 
the repeated interactions of section [2T5l as difference equations. Clearly we have 

AU{1) = U{1) - U{1 - 1) = {Ml - I) U{1 - 1), C/(0) = /. 

We will need to convert this expression into a more convenient form, however, in order 
for it to be of use in calculations. 

Recall that A/;, and hence also Mi — I, is only a function of the initial system and 
of the ^th time slice of the field; indeed, this is easily read off from Eq. (j2.7p . But any 
operator X G ^ iS" ^ can be written in a unique way as 

X = X^(g) cr+cr_ +X+ ®C7++X- (^a^+X° e 

as cr+fT-, (T+, (T_, / is a linearly independent basis for By the same reasoning we 
find that we can write without loss of generality 

Mi-I = M^ AA(0 + M+ AA*{1) + M' AA{1) + M° At{l), 

so that M°'+^^^* e ^0 are uniquely determined by Li_2.3 and A. The functional 
dependence of M°'+^^^* on ii.2,3 and A is complicated and we will not attempt to give 
a general expression; in the examples which we consider later it will be straightforward 
to calculate these matrices explicitly. 

Remark 4.4. It should be emphasized that whereas ^1,2,3 G do not depend 
on the time step size A, M°'+'~'^ are in fact functions of A. As usual, we have chosen 
the scaling of the various parameters with A so that the limit as A ^ gives meaningful 
results; in this case, it can be shown that the matrix elements of M°''^'~'^ converge 
to finite values as A ^ 0, and the limiting matrices can be expressed explicitly in 
terms of Li^2.3 (see 02] and other references in section [TOl) . 

We now obtain the following quantum stochastic difference equation: 

AU{1) = {M± AA(0 + M+ AA*{1) + AA{1) + M° At{l)} U{1 - 1), (4.1) 

with U{0) = I. It is this form of the equation that will be the most useful to us, as 
we can apply stochastic calculus techniques to manipulate it. We could equivalently 
express it in terms of discrete quantum stochastic integrals: 

U{1) =1+ {M^U-) ■ k{l) + {M+U-) ■ A*{1) + {M-U-) ■ A{1) + {M°U-) ■ t(l), 

where U-{1) = U{1 — 1). This is the way that (quantum) stochastic differential 
equations are defined in continuous time. 

To get some familiarity with calculations with the quantum Ito rule and Lemma 
14.31 let us calculate the time evolution of the expectation V{ii{X)) for an atomic 
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observable X g By linearity of the state we cleariy have AP(j;(X)) = P(Aj/(X)), 
so we are going to calculate Aj;(X) — A{U{1)* XU{1)) using the quantum Ito rule 
and then calculate the expectation of the resulting expression. First, note that 

AU{1)* = U{1 - 1)* {M^* AA(0 + M+* AA{1) + M-* AA*{1) + M°* Atil)} . 

Then we calculate using the quantum Ito rule 

A{U{iyxU{l)) -(•••) AA(/) + {■■■) AA*{1) + {■■■) AA{1) 
+ U{1 - 1)* {M+*XM+ + M°*XM° + M°*X + XM°) U{1 - 1) At(l). (4.2) 

We didn't bother to calculate the AA, AA, AA* terms: by Lemma 14.31 these vanish 
anyway when we take the expectation with respect to P, so we only need the At term. 
Hence we obtain for any X G SSq 

^m.n.-(£m)). («) 

where the discrete Lindblad generator £(•) is defined by 

C{X) = M+*XM+ + \^ M°*XM° + M°*X + XM°. (4.4) 



Remark 4.5. Let Xt be a classical, continuous time Markov diffusion, and 
define the semigroup Pt{f) = Ef{Xt). Then (under sufficient regularity conditions) 
dPt{f)/dt — Pt{Cf), where £ is the infinitesimal generator of the Markov diffusion. 
Eq. (|4.3p is strongly reminiscent of this formula, and indeed the continuous time 
version of the formula plays an equivalent role in quantum probability. In fact, the 
semigroup property (expressed in "differential" form (|4.3p ) suggests that we should 
interpret our repeated interaction model as a "quantum Markov process" , which can 
be given a precise meaning but we will not do so here. 

Note that we are not allowed to choose M^'+'~'° arbitrarily: these matrices 
must be chosen in such a way that the solution U (l) of the corresponding difference 
equation is unitary. Obtaining M^'+'^'° directly from Li,2.3, as we will do in the 
examples below, ensures that this is the case. One could establish general conditions 
on M^'+'~'° to ensure that U{1) is unitary, but we will not need these in the following. 
There is one particular necessary condition, however, that we will often use. 

Lemma 4.6. For any repeated interaction model C{I) = 0. 

Proof. As U{1) is unitary, = U{l)*U{l) = I for any /. By Eq. (g^]) we have 

But this must hold for any initial state p, so C{I) =0 identically. □ 

4.3. Examples. Let us now return to the examples of section 12.61 We will 
calculate the difference equations for the unitary evolution U{1) explicitly. 

Spontaneous emission. Recall that in this case (with 2k, = I) 



Ml = exp (ct_ AA*{1) - a+ AAil)) 
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As Ml is of the form P O 7®'"! Q O 7®'=-', we can read off the form of M; by 
calculating the matrix exponential of 

/O 0\ 

A 

-A 

\0 0/ 

Performing this calculation explicitly, we get 



/I 













cos A 


sin A 








— sin A 


cos A 













1/ 



\ , I (8) fT+cr_ + sin A cr_ (g) cr+ — sin A cr+ (g) cr_ 

cos Ay ^ ^ ^ 

= (1 — cos A) Gz ® a+ff- + sin A cr_ (g) ct+ — sin A C7+ (g) a_ + 



cos A 0\ 

iJ®^-'^+ 

cos A 0^ 
1> 



where we have used <j-<j+ = I — a+a-. Hence we obtain 

Mi = {l- cos A) az AA(0 + ^ (a_ AA* (/) - ct+ AA(Z)) + ~ a+(T_ + 1. 

A A 

We can now immediately read off the coefficients in the quantum stochastic difference 

equation for the spontaneous emission model: 



= (1 - cosA)f7z, 



, ^ , sin A 



M- 



sin A 



M° 



A 

cos A — 1 



(4.5) 



A " A2 

Dispersive interaction. Recall that in this case (with g = 1) 

Ml = exp(cr^(A^*(/) - A^(0)) . 

We proceed as in the case of spontaneous emission. Starting from 



A (cr+ — cr_ ) 



( " 





A 


\ 











-A 


-A 











^0 


A 





0^ 



we calculate the matrix exponential 
/ cos A 

pTz^X (<7+— 0-_) _ 



Hence we obtain 

Ml 





— sin A 

V 

sin A 






sin A 




cos A 





— sin A 





cos A 





sin A 





cos A j 



— cosA7 + sinACT2 (ct_|_ — cr_). 



a,(AA*(0-AA(0) + 



cos A — 1 



A ■ A2 

We can now immediately read off the coefficients in the quantum stochastic difference 
equation for the dispersive interaction model: 

sin A 



M^ 



0, 



M+ 



-M- 



M° = '^^I. (4.6) 
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Fig. 4.1. The expectation of the energy of a spontaneously emitting atom as a function of time. 
The initial state p was chosen in such a way that the atomic energy takes its maximal value +hLDo/2 
with unit probability at t = 0, i.e. p{x) = Tr[(T+(T„X]. The expected energy decays exponentially, 
until it reaches its minimum value — fkiJo/2. The time scale used for the calculation is = 300. 



A first simulation. To get some idea for the mean behavior of our two examples, 
let us calculate the expectation of the energy observable F{ji{H)) as a function of time. 
The expectation of any atomic observable is given by Eq. (|4.3|) . To simulate Eq. (|4.3p 
directly, we can use a standard technique which will be used several times in the 
following. Note that P{ji{X)) is a linear function of X. Hence we can always find a 
2x2 matrix t; such that Tr[TiX] = F{ji{X)) for every X. Substituting this expression 
into Eq. ()4.3|) . we obtain explicitly a recursive equation for Tf. 

where tq is the density matrix corresponding to the initial state p: p{X) = Ti[tqX] for 
every X. This equation is called the (discrete) master equation, and plays a similar 
role to the forward Kolmogorov (or Fokker-Planck) equation in classical diffusion 
theory. The recursion is easily implemented in a computer program, and the solution 
allows us to calculate F{ji{X)) for any X. For the spontaneous emission case, the 
expected energy P{ji{H)) — Tr [rjiJ] is plotted in Figure [57TJ 

We have used this opportunity as an excuse to introduce the master equation; in 
this simple case, however, we can obtain the result of Figure HTTl much more directly. 
Let us calculate C{H) = {hujQ/2) £{az) in the case of spontaneous emission: 

2 sin^ A sin^ A , ^, 
^'^z) = j^a+a- = — (cr^ + /). 

Hence we have 

nM'^z)) + 1 = AP(jj(a,)) +P(j,_i(a,)) + 1 = cos2 A(P(j,_i(a,)) + 1). 
Recursing this relation, we find that the expected energy decays geometrically: 

njim = ~^ + [piH) + ^) (cos^Ay. 

This expression coincides exactly with the plot in Figure W7l\ 

In the dispersive case, the time dependence of the mean energy is trivial: it is not 
difficult to check that in this case C{H) = 0. This is not too surprising, as we have 
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suppressed the atom's ability to radiate its excess energy by placing it in an optical 
cavity. In fact, using A// = exp(cr2(AA*(/) — AA{1))) we find that 

ji{H) = Ml ■ ■ ■ Mi%M;HMiMi^i ■■■Mi=H. 

Evidently not only the expected energy, but even the energy as a random variable, is 
conserved in the dispersive case. 

5. The martingale method. In this section we will provide a solution to the 
filtering problem using martingale methods. The key tools we need are minor varia- 
tions on two classical theorems. First, there is a theorem of Doob which states that 
any adapted process can be decomposed as the sum of a predictable process and a 
martingale; we extend this theorem slightly to a class of nondemolition processes. Ap- 
plying this theorem to the conditional expectation 7rz(X), we can immediately identify 
the corresponding predictable part. Next, the martingale representation theorem tells 
us that any martingale can be written as a stochastic integral of some predictable pro- 
cess. Consequently we can identify also the martingale part of iti{X) in terms of a 
stochastic integral. The recursive nature of the solution is the motivation behind this 
procedure: the predictable parts will exactly turn out to be functions of the filter in 
the previous time step, whereas the increment in the discrete stochastic integral will 
be directly related to the observation increment. 

For excellent introductions to the use of martingale methods in classical filtering 
theory we refer to [29l [60] ; definitive treatments can be found in [65l |56] . Martingale 
methods for quantum filtering were introduced in [1^ and used in [T8l [20] , The 
treatment below is a discretized version of the latter. 

5.1. The classic theorems. Recall that {^} is the commutative filtration gen- 
erated by the observations. 

Definition 5.1. A quantum process X is called nondemolished (by the observa- 
tions Y) if X{1) is in the commutant of ^Vi for < I < k. A nondemolished quantum 
process H is called an nd- martingale (with respect to the filtration '3^o<i<k) if 

V{H{l)\?Kn) = P(iJ(m)|^„) VO < m < / < k. 

A quantum process A is called -predictable if A(V) is in '3^i-\ for 1 < I < k, and A is 
called ^^-adapted if A{1) is in 'Wi for 1 < I < k. An nd-martingale that is additionally 
W -adapted is called a ^^-martingale or simply a martingale. 

The nondemolition requirement ensures that the conditional expectations in the 
definition of an nd-martingale are well defined: C implies D so 

H{1) g implies H{1) e Note also that by construction, any ^^^-predictable or 
^ -adapted process is a function of the observations and is hence commutative. 

Remark 5.2. The concept of a ^^-martingale, being ^^^-adapted and hence a 
classical process, coincides with the classical notion of a martingale. An nd-martingale 
("nd" for nondemolition) is used in a slightly broader sense; we will encounter an 
example below of a commutative nd-martingale that is not a martingale. It will be 
convenient to use this terminology, as we will use the following theorem. 

Theorem 5.3 (Doob decomposition). Let X be a quantum process nonde- 
molished by the observations Y . Then X has the following Doob decomposition 

X{1) = X{0) + A{1) + H{1), l^l,...,k 

where H is an nd-martingale (w.r.t. '3^{)<i<k) null at 0, and A is a "3^ -predictable 
process null at 0. Moreover, the decomposition is unique modulo indistinguishability. 
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The proof is very similar to its classical counterpart [55]. 

Proof. Suppose X has a Doob decomposition as in the Theorem. Then 

p(x(0-x(/-i)|^^_i) 

= F{H{1) - H{1 ~ iWi-i) + P(A(0 - A{1 - iWi-i) = ^(0 - A{1 - 1), (5.1) 
where we have used that H is an nd-martingale and that A is predictable. Hence 

^P(X(m)-X(m-l)|^^„_i), A(0) = 0. (5.2) 

m— 1 

For any nondemolishcd process X , define the predictable process A{1) as in Eq. (|5.2p . 
and define H{1) = X{1) — X{0) — A{1). Then it is easily verified that H is an nd- 
martingale, hence we have explicitly constructed a Doob decomposition. 

To prove uniqueness, suppose that X{1) — X{0) + A{1) + H{1) is another Doob 
decomposition. Then A{1) — A{1) — H{1) — H{1) for 1 <l <k, and hence 

A{l)-A{l)^nA{l)-A{lWi-,) 

= P{H{1) - H{lWi-i) = nH{l - 1) - H{1 - iWi-i) 

^ ¥{A{1 - 1) - A{1 - l)|^-i) - A{1 - 1) - A{1 - 1) (5.3) 

where we have used predictability of A, A and the nd-martingale property of H, H. 
But as A{0) = A{0) = 0, we obtain by induction ^ A{1) - A{1) = H{1) - H{1) 
for 1 < I < k. Hence A = A and H = H with probability one (as the conditional 
expectations in Eq. (|5.3p are only defined up to a choice of version). □ 

Remark 5.4. It should be noted that the Doob decomposition depends crucially 
on the choice of filtration ^ which is demonstrated by the following trivial example. 
Consider the (commutative) filtration {J'l} with .^i — alg{/} for any I. As contains 
only multiples of the identity, the commutant is the entire algebra ^' — ./# ® Wk- 
Hence any process X is nondemolishcd by J^, and as P(X(/)| J^„i) = IP(-'^(0) foi' ^"^y 
l,m any process with constant expectation P{X{1)) = ¥{X{m)) is an nd-martingale 
with respect to {J^i} (but not necessarily an ^-martingale!) Using as the filtration, 
we obtain the Doob decomposition A{1) = P{X{1) - X{0)) and H{1) X{1) - X{0) - 
F{X{1) — X{0)) for any process X. Clearly this decomposition is different than the 
Doob decomposition with respect to ; but note that H{1) is not an nd-martingale 
with respect to so uniqueness is not violated. The moral of the story is that we 

have to be careful to specify with respect to which filtration we decompose a process. 
In the following, this will always be the filtration {^^} generated by the observations. 

Let X e Applying the Doob decomposition to tti{X) and Y{1) gives 

7Ti{X)=p{X) + B{l) + H{l), i.e. Atti{X)^ AB{l) + AH{l), 
Y{l)=Y{0) + C{l) + Y{l), i.e. AY{l)^AC{l) + AY{l), ^^'"^^ 

where B and C are predictable processes null at and Y and H are -martingales null 
at 0. The process Y is called the innovating martingale. In the next two subsections 
we will investigate the processes C and Y in more detail for both the counting and 
homodyne detection cases. 

Lemma 5.5. The predictable process B in the decomposition ofTTi(X) is given by 

AB{1) ^ Tri^i{C{X)) At{l) l<l<k, 
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where C is the discrete Lindblad generator of Eq. (14. 4p . 

Proof. By Eq. jSl]), we have AB(l) = P(A7r;(X)|^_i) = P(Ajz(X)|^_i). To 
calculate the latter, let K be an element in Using Eq. (|4.2p . we obtain 

P(i^Aj,(X)) =P(if j,_i(/:(X))At(0) =P(/^7r,_i(/:(X))A<(0). 

As this holds for any K E "S^i-i, and as 7r;_i(£(X)) e the statement of the 

Lemma follows from the definition of the conditional expectation. □ 
At this point, we have the expression 

AniiX) = Tri^i{C{X)) At{l) + AH{1). 

This is almost a recursive equation: what we would like to do is write something like 
AH{1) = /(tt;-! (•)) Ay(?), as in that case we could use this equation to calculate 
TTi{-) using only and AY{1). The problem is that /(7r;_i(-)) AF(/) does not 

define a martingale; but /(7r;_i(-)) AF(/) does! This suggests that we should try to 
represent i? as a discrete stochastic integral with respect to the innovating martingale 
Y{1). The martingale representation theorem shows that this is always possible. 

Theorem 5.6 (Martingale representation). LetY be the innovating martin- 
gale and let H be a -martingale null at 0. Then there exists a -'predictable process 
S such that AH{1) = AY(l), I ^ 1, . . . , k modulo indistinguishability. 

The following proof is reminiscent of [82j pp. 154-155], but the details of the 
argument are a little more delicate in our case. 

Proof. As all the observables in the theorem are contained in the full observation 
algebra which is commutative, this is essentially a classical problem. It will be 
convenient for the proof to treat it as such, i.e., applying the spectral theorem to 
(^%,P) gives the classical probability space {fl,!F,P), the filtration {Wi} gives rise 
to the classical filtration {yi}, and we will write yi = i{Y{l)), jji = l{Y{1)), and 
hi = l{H{1)). It will be convenient to write ivi = i,{AY{l)) = yi — yi-i- 

We will make fundamental use of the following fact: uji takes one of two values 
{lu+,lu-} for every I = l,...,k. To see this, recall that AY{1) = U (1)* AZ {l)U (1) 
where AZ{1) is one of AA{1) + AA*{1) or AA{1) (in fact, any observable AZ{1) of 
the form {Z)i G Wk, Z £ ^ would do, provided Z is not a multiple of the identity.) 
Hence AZ{1) has two distinct eigenvalues, and as unitary rotation leaves the spectrum 
of an operator invariant so does Ay(^). It follows that l{AY{1)) is a two-state random 
variable. We will write pf = P(w/ = uj^\yi-i) for the conditional probability that 
uji — uj^ given that we have observed oji, . . . , 

Now recall that ^ is the algebra generated by AY{i), i — 1, . . . ,1. Hence every 
3^/ -measurable random variable can be written as a function of w;, i — I, . . . ,1. In 
particular, yi = j/i(wi, . . . , w/) and hi = /i;(wi, . . . , w/). We would like to find a 
^;(a;i, . . . ,cj;_i) (independent of uji, hence predictable) that satisfies 

Ahi{ui, ...,uji)^ 6(t^i, • • ■ ,^i-i) Ayi{uji, . . .,uji). (5.5) 

To proceed, we split fl into three disjoint subsets ^1.2,3 and define the random variable 
S.I separately on each set. 

Case 1: 17i = {w € il : Ayi{uji, . . . ,uji-i,uj'^) ^ 0}. Let us suppose that ^/ 
exists. Then evidently on Qi 



6(wi, • . .,UJl-l) 



Ahijuji, ...,uji) 
Ayi{uji, ...,uji)' 
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Existence of is thus verified by construction if we can show that the right hand side 
is independent of uji. To this end, we express the martingale property of hi as 

Ahi{uJi, . . . ,a;/_i,Lj+)p,+ (LJi, . . . ,cj;_i) + Ahi{u!i, . . . ,uji^i,uj^)pf {uji, . . . ,cj;_i) = 

a.s., where the left hand side is simply the expression for Ep(A/i;|3^;_i). Similarly, 

Ayi{LL>i, . . . . . . + Ayi(wi, . . . ,ui-i,uj^) p,^ [uJi, . . . ,uji-i) = 

a.s., as jji is a martingale. Note that necessarily pf^ ^ a.s. on fli. Hence we obtain 

A/l;(wi, . . . ,a;;_i,LJ+) 

Ayi{uji, . . . , Wi-i, w+) 

Dividing the first by the second expression, the independence of ^; from uji follows. 

Case 2: = {w G : p'^ € {0, 1}}. Using the martingale property of hi and yi 
as above, we conclude that on 1^2 we have Ahi = Ayi = a.s. Hence Eq. (|5.5|) holds 
regardless of the value we assign to ^i . 

Case 3: ft^ = U ft2)- We will show that P(fi3) = 0, so that we do not 

need to worry about defining ^i on this set. As on fl^ we have pf^ ^ but one 
of Ayi{LUi, . . . ,uJi^i,u}^) = 0, using the martingale property as above allows us to 
conclude that Ay; — a.s. on fia. Recall that uji = Ayi = Aci + Ayi where Aci 
is predictable. Then uji ~ Aq(wi, • • • a.s. on Q^. But this would imply that 

uJi = Ep((jj;|3^;_i) = lli^pI^ +Lj~p'j^ , and as pf^ ^ we would be able to conclude that 
uJi ^ uj^ ■ Hence we have a contradiction. 

We have now shown how to define ^i that satisfies Eq. (|5.5p except possibly on a 
set of measure zero. Setting S(Z) — the theorem is proved. □ 

Though the proof of the discrete martingale representation theorem is in princi- 
ple constructive, it is not advisable to follow this complicated procedure in order to 
calculate the predictable process S. Instead we will calculate S using a standard trick 
of filtering theory, and it will turn out to depend only on the conditional expectations 
in the previous time step. Putting everything together, we obtain a recursive relation 
with which we can update our conditional expectations of atomic operators given the 
conditional expectations at the previous time step and the observation result at the 
present time step. This recursion is called the discrete quantum filtering equation. 
As the predictable processes C and S depend on the nature of the observations, we 
consider separately the homodyne and photon counting cases. 

5.2. Homodyne detection. Let us first consider a homodyne detection setup, 
i.e. an experimental setup that allows us to observe 

Y^{1) = U{l)*{A{l) + A*{l))U{l), 0<l<k. 

We begin by finding the predictable process C in the Doob decomposition of . 
Lemma 5.7. The predictable process C in the decomposition ofY-^ is given by 

AC{1) = TTi^i{M+ + M+* + \^ M°* M+ + \^ M+* M°) At{l), l<l<k. 



= -Ahi{uji, . . .,L^Ji^i,uj ) 



-Ayii^uji, 



. . . 

(CJI, . . .,UJl-l) 
pt{uJi, ■ ■ . ,t^i-l) 
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Proof. By Eq. dUI]), we have AC{1) = P(Ar-^(/)|^_i). To calculate the latter, 
let K be an element in we would like to find an expression for ¥{KAY^(l)). 

To this end, we calculate using the discrete quantum Ito rules 

A{{A{1) + A*{l))U{l)) = (•••) AA(/) + (...) AA{1) 

+ {!+ iA{l -1) + A*{1- 1)) M+ + \^M°} U{1 ~ 1) AA*{1) 

+ {iA{l - 1)+A*{1 - 1))M° +M+} U{l-l)At{l), (5.6) 

where we have retained the relevant terms. Consequently, we calculate 

A{U{l)*{A{l) + A*{l))U{l)) = (•••) AA(0 + {■■■) AA*{1) + {■■■) AA{1) 
+ U{1 - l)*L{I){A{l -l) + A*{l- l))U{l - 1) At{l) 

+ ji-i{M+ + A/+* + X^M°*M+ + \^M+*M°) At{l), 

where C{X) is the discrete Lindblad generator. But £(/) vanishes by Lemma 14.61 
Hence we obtain 

fiKAY^il)) = ¥{K ji_i(A/+ + M+* + \^M°*M+ + A^A/+*Af°) At{l)) 

using Lemma 1431 or equivalently (using K e '3^i-i) 

F{KAY^{1)) = P(i^ 7r;„i(Af + + M+* + \^M°*M+ + \^M+*M°) At{l)). 

As this holds for any K e ^^-i, and as tti-i{- • • ) G ^^-i, the statement of the Lemma 
follows from the definition of the conditional expectation. □ 

From Theorem 15. 6l we know that AH = EAY for some predictable process 2. It 
remains to determine S; we approach this problem using a standard technique. Since 
¥{ji{X)Y{l)\i^i) = Tri{X)Y{l), the uniqueness of the Doob decomposition ensures 
that ji{X)Y{l) and tti{X)Y{1) have equal predictable parts. We write D and E for 
the predictable processes in the Doob decomposition of ji{X)Y{l) and tti{X)Y{1), 
respectively. Solving the equation AD = AE will then allow us to determine S. 

Lemma 5.8. For any X e Mo, define 

J{X) = XM+ + M+*X + \^M°*XM+ + \^M+*XM°, 

so that AC{1) — '7r;_i (j/(/)) At{l). Then we can write 

^ 7:i-i{J{X)) - 7ri_i(X + A^£(X)) ni^i{J{I)) 
' I - X'^i-iiJil)? 

Proof. We begin by evaluating AD. For any K E ^-i we want to calculate 
P{A{ji{X)Y{l))K) ^ f>{A{U{l)*X{A{l) + A*{l))U{l))K). We proceed exactly as in 
the proof of Lemma 1 5. 71 Using Eq. (15. 6p and the quantum Ito rules, we obtain 

A{Uil)*X{A{l) + A*{l))U{l)) = (•••) AA(/) + {■■■) AA*{1) + {■■■) AA{1) 
+ U{1 - \)*C{X){A{1 -l) + A*{l- l))U{l - 1) At{l) 

+ ji-i{XM+ + M+*X + X^M°*XM+ + X^M+*XM°) At{l). 

It follows in the usual way that 

AD{1) ^ 7n^i{C{X))Y{l ~ I) At{l) 

+ ■Ki-i{XM+ + M+*X + \^M°*XM+ + \^M+*XM°) At{l). 
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We now turn our attention to A£'. First note that the Ito rule gives 

A{7Ti{X)Y{l)) - {Atti{X)) Y{1 - 1) + 7Ti_,{X) AY{1) + Atti{X) AY{1). 

By uniqueness of the Doob decomposition, AE is the sum of the predictable parts 
^£'1,2,3 of the three terms on the right hand side. Let us investigate each of these 
individually. The first term can be written as 

{Atti{X)) Y{1 - 1) = 7ri_,{C{X)) Y{1 - 1) At{l) + Y{1 - 1) AH{1). 

It is easily verified, however, that Y{1 — 1) AH{1) inherits the martingale property 
from H. Hence by uniqueness of the Doob decomposition, we obtain AEi{l) = 
7r;_i(£(X)) Y{1 — 1) At{l). Moving on to the second term, we write 

7Ti-i{X) AY{1) = 7r,_i(X) AC{1) + n^iiX) AY{1). 

Similarly as above we find that 7r/_i(X) Ay(/) inherits the martingale property 
from Y, so that evidently Ai?2(0 = AC{1) (where AC{1) is given explicitly 

in Lemma [S77|) . It remains to deal with the third term. To this end, let us write 

AiTiiX) AY{1) ^ {ni^,{C{X)) At{l) + E{1) AY{1)} {AC{1) + AY{1)). 

As before, processes of the form ^;X(/)Ay'(/) with ^^-predictable X inherit the 
martingale property from Y . Thus we only need to retain the predictable terms: 

AniiX) AY{1) = (•••) AY{1) + tti-i{£{X)) At{l) AC{1) + {AY{l)f. 

Similarly, we expand {AY{1))'^ as 

{AY{iyf = (...) AY{1) + {AY{l)r - (AC(/))2, 

where we have used AY{1) — AC{1) + AY{1). But using the Ito rules we calculate 
(Ay(/))2 U{l)*{AA{l) + AA*{l))^U{l) = At{l). Hence we can read off 

AEsil) = TTi-i{£{X)) At{l) AC{1) + E{1) {At{l) ~ {AC{1))^). 

But recall that ji{X)Y{l) and tti{X)Y{1) have equal predictable parts. Hence setting 
AD{1) = AEi{l) + AE2{1) + AE:i{l) and solving for the Lemma follows. □ 
Putting everything together we get the following discrete quantum filtering 
equation for homodyne detection 



Am{X) ^ ^i^^{C{X)) At{l) 

ni-i{J{X)) ~ iri^ijX + X^CjX)) m-i{J{I)) 
/-A27r,_i(J(/))2 



(Ay(0-^i-i(^(/))At(0). 



5.3. Photodetection. We now turn our attention to a setup where we are 
counting photons in the field, i.e. we are observing 

Y^{1) ^U{l)*A{l)U{l), 0<l<k. 

The procedure here is much the same as in the homodyne detection case. 

Lemma 5.9. The predictable process C in the decomposition ofY^ is given by 

AC{1) ^ 7r,_i (M+*M+) At{l), l<l<k. 
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Proof. By Eq. dm]), we have AC(/) = P(Ay^(0|^_i). To calculate the latter, 
let K be an element in we would like find an expression for P{KAY-^{1)). To 

this end, we calculate using the discrete quantum Ito rules 

A(A(OC/(0) = (---)AA(0 + (---)AA(0 

+ (/ + A(/ - l))M+U{l - 1) AA*{1) + A{1 - l)M°U{l - 1) At(/), (5.7) 

where we have retained the relevant terms. Consequently, we calculate 

A(f/(0*A(/)[/(0) = (•••) AA(0 + {■■■) AA*{1) + {■■■) AA{1) 

+ U{1 - l)*C{I)K{l - l)U{l - 1) + j,_i(Af+*M+) At(0, 

where C{X) is the discrete Lindblad generator. But C{I) vanishes, and the Lemma 
follows by the usual argument. □ 

Next, we determine the predictable process such that AH{1) — E{1) AY{1). 

Lemma 5.10. The process is given by 

E{l) = {I-X\i.iiM+*M+))-^ 

Proof. We begin by finding the predictable process AD in the Doob decompo- 
sition of ji{X)Y{l). For any K e ^-i we want to calculate V{A{ji{X)Y{l))K) = 
P{A{U{l)*XA{l)U{l))K). Using Eq. ^Ji) and the quantum Ito rules, we obtain 

A(f/(0*XA(/)[/(0) = (•••) AA(/) + {■■■) AA*{1) + {■■■) AA{1) 

+ U{1 - 1)*£{X)A{1 - 1)U{1 - 1) At{l) + ji_i{M+*XM+) At{l). 

It follows in the usual way that 

AD{1) = 7ri^i{C{X)) Y{1 - 1) At{l) + Tri-i{M+*XM+) At{l). 

We now turn our attention to the predictable process AE in the Doob decomposition 
of TTi{X)Yi. First note that the Ito rule gives 

A{7n{X)Y{l)) = (ATTiiX)) Y(l - 1) + n-i{X) AY{1) + ATn{X) AY{1). 

By uniqueness of the Doob decomposition, AE is the sum of the predictable parts 
A£'i^2.3 of the three terms on the right hand side. As in the proof of Lemma [5?8l we 
find that AEi{l) = 'Ki-i{C{X))Y {I - I) At{l) and AE2{1) = 7r,_i(X) AC(/) (where 
AC {I) is given explicitly in Lemma |5.9[) . To deal with the third term, we write 

A^i{X) AY{1) = {iTi^^{C{X)) At{l) + S(0 Ay(/)} {AC{1) + AY{1)). 

As before, processes of the form 'Y^iX{l) AY{1) with ^^-predictable X inherit the 
martingale property from Y . Thus we only need to retain the predictable terms: 

A^i{X) AY{1) = (•••) AY{1) + ^i-i{C{X)) At{l) AC{1) + E{1) {AY{1))\ 

Similarly, we expand {AY{1))'^ as 

(Ay(/))2 ^ (. . . ) AY{1) + {AY{l)r - (AC(0)^ 



(M+*A:Af+) 



ni-i{M+*M+) 



TTl^liX + X''C{X)) 
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where we have used AY(l) — AC{1) + AY{1). But using the Ito rules we calculate 
(Ar(/))2 = U{l)*AMlfU{l) = AY{1). Hence 

{AY{l)f = (...) Af (0 + AC{1) - {AC{1))\ 

and we can read off 

AE:i{l) = Tii_i{C{X)) At{l) AC{1) + {I - AC(/))AC(0. 

But recall that ji{X)Y{l) and tti{X)Y{1) have equal predictable parts. Hence setting 
AD{1) = AEi{l) + AE2{1) + AE-iil) and solving for E{1), the Lemma follows. □ 
Putting everything together we obtain the following discrete quantum filtering 
equation for photon counting 

AniX) = Tri^i{C{X)) At{l) + (/ - \\i^i{M+*M+))-'x 
'■Ki-i{M+*XM+) 



TTl^l{M+*M+) 



-Tii^^{X + \'C{X)) 



{AY{1) - 7r,_i(M+*M+) At{l)). 



5.4. How to use the filtering equations. The filtering equations of sections 
and 15.31 mav seem a little abstract at this point; 7rj(X) is some observable in 
the algebra ^ (8) #fc, and it appears that we would need to know 7ri_i(Z) for every 
Z e SSq^ in addition to the observation increment AY {I), in order to be able to 
calculate t^i{X) for arbitrary X. The equations are much less abstract than they 
might seem, however. First of all, recall that both 'rTi{X) and AY{1) are elements of 
the (commutative) observation algebra in fact, the filtering equations live entirely 
within this algebra. Hence these are just classical equations in disguise (as they should 
be!); we could write explicitly, e.g. in the homodyne detection case. 



Ai{Tn{X)) = L{n-i{C{X))) a2 + {i{AY{l)) - t(^,_i(j7(/))) A^)x 

L{n-i{J{X))) - i{m-i{X + \^C{X))) i{7ri^,iJ{I))) 



(5.8) 



using the *-isomorphism l obtained by applying the spectral theorem to (^%,P). For 
any < I < k and X G l{7Ti{X)) is a random variable that is a function of the 
random process i,{AY{i)) up to and including time I. But an elementary property of 
the conditional expectation is that ni{X) is linear in X; hence l(t:i{X)) is also linear 
in X. This means we can always write L{Tri{X)) = Tt[piX], where pi is a (random) 
2 X 2-matrix (as ./# is two-dimensional). We obtain the following recursion for pf. 

w ^ ^2 , Jjpi-i) - Tr[J(p;-i)] {pi-i + X'C{pi-i)) ^, 

+ l-A^Trm,._0]^ _ 

(.(Ay(/))-Tr[^(p^_i)]A^), (5.9) 

where we have written 

C{p) ^ M+pM+* + M°pM°* + M°p + pM°\ 
J{p) = M+p + pM+* + \^M+pM°* + X'^M°pM+*. 

Po is simply the density matrix corresponding to the initial state p{-), i.e. p{X) = 
Tr[po^] for every X G ./#. The matrix pi is called the conditional density matrix and 
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contains all the information needed to calculate tt; (X) for every X g Furthermore, 
Eq. (|5.9p is a simple nonlinear recursion for 2 x 2-matrices. At any time step / we only 
need to remember the 2 x 2-matrix pi; when the {I + l)th observation l{AY{1 + 1)) 
becomes available, which takes one of the values ±A, we simply plug this value into 
Eq. (|5.9p and obtain the updated matrix pi+i. Such a recursion is very efficient and 
would, if necessary, be easily implemented on a digital signal processor. 

As filtering equations are entirely classical, there is no real need to make the 
explicit distinction between their representation in terms of classical random variables 
on a probability space vs. elements of the observation algebra. We urge the reader 
to always think of sets of commuting observables as random variables: this is implied 
by the spectral theorem, and is at the heart of quantum mechanics! Eq. (|5.8p is 
notationally tedious and completely unnecessary, as it does not add anything to the 
filtering equation as we have already written it in section 15.21 In some cases, e.g. in 
the proof of the martingale representation theorem, it is convenient to use explicitly 
the structure of the underlying probability space; but in much of this article we will 
not make an explicit distinction between random variables and observables. 

A similar story holds for the photodetection case; we leave it up to the reader to 
calculate the associated recursion for the conditional density matrix. 

5.5. The Markov property and Monte Carlo simulations. The filtering 
equations that wc have developed take as input the observation process obtained from 
the system. Though this is precisely how it should be, one would think that further 
investigation of the filters can not proceed without the availability of typical sample 
paths of the observations from some other source, be it an actual physical system 
or a direct computer simulation of the underlying repeated interaction model. It 
is thus somewhat surprising that we can actually simulate such sample paths using 
the filtering equation only, without any auxiliary input. This is due to the Markov 
property of the filter, which we will demonstrate shortly. We can use this property 
of the filtering equations to perform Monte Carlo simulations of both the sample 
paths of the observation process and sample paths of the filter itself (called "quantum 
trajectories" in the physics literature). In addition, the Markov property is key for 
the development of feedback controls, as we will see in sections [7]-[9l 

We will consider below the homodyne detection case, but the photodetection case 
proceeds identically. Set Ay; — l{AY{1)). Recall that the homodyne detection signal 
Ay; takes one of two values ±A for every I. Suppose we have observed y^ up to and 
including time I — I', we would like to be able to calculate the probability distribution 
of Ay; using this information. This calculation is carried out in the following lemma. 

Lemma 5.11. We have 

P[Ay; = ±A I yi^,] = p{Ayi = ±A; p;_i), (5.10) 

where 

p(Ay-±A;p)-i±^Tr[J(p)] (5.11) 

depends only on the filter in the previous time step. 

Proof. Let = P[Ay; = ±A | 3^;_i] be the probability that the observation in 
the next time step Ay; takes the value +A. Using the Doob decomposition we have 



AY{1) ^ AC{1) + AY{1) ^ Tri-i{J{I)) At{l) + AY{1), 
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where Ay(/) is a martingale increment, so that 

Xp+ - A(l - p+) = /,(P(Ar(0|^-i)) = L{ni-i{J{I))) \\ 

Thus 

Pt = \ + \ Trh_i^(/)] = i + ^ Tr[J(p,_i)] 

depends only on the filter in the previous time step. □ 

With this distribution in hand, we can also calculate the statistics of the filter in 
the next time step, and we can prove the Markov property by recursing this procedure. 

Lemma 5.12. The filter pi satisfies the Markov property: 

Ep(5(ft)|a{po, . . ■,Pi}) = Ep{g{p,)\a{pi}) Wl<j<k. 

Proof. Using the recursion ()5.9p . pj can be written as a deterministic function 
f{pj-i, Ayj) of pj-i and Ayj. By the martingale property of the innovation process 

P(Ay, - ±X\a{yi, . . . = ^ ± ^ Tr[J(p,_i)], 

which is only a function of Pj-i- As a{po, . . . , Pj-i} C cr{yi, • . • , J/j-i}, we obtain 

P{Ay, = ±X\<j{po, . . . = ^ ± ^ Tr[J(p,_i)]- 

Hence for any function g 



^p{9iPj)W{Po, ■ ■ ■,Pj'i}) = Ep(5'(pj_i, Ayj)|cr{po, . . ■,Pj-i}) 



J2 9'iPj-uaX) i+a^Tr[^(p,_i)] 



ae{~lS} 

where we have written g'[p,w) — g{f{p,w)) for the function g composed with the 
one-step filter recursion. As the right hand side is a function of pj-i only, we have 

Ep(5(Pj)k{Po, • • ■ = Ep(Ep(g(pj)|cr{po, • • ■ , Pj-i}) I cr{pj_i}), 

from which we conclude that 

Ep(5(Pj)k{Po, • ■ ■ , Pj-i}) ^ ^p{giPj)W{pj-i}). 

But setting 'Eip{g{pj)\a{po, . . . = /i(pj-i), we can repeat the argument giving 

Ep(/l(pj_i)|cr{po, • ■ • , Pj-2}) = Ep{h{pj^i)\a{pj^2})- 

From the definition of h{-) we immediately obtain 

Ep(.9(pj)|cr{po, • ■ ■,Pj-2}) = Ep{g{pj)\a{pj^2})- 

Recursing the argument gives the Markov property. □ 

It is now straightforward to turn this procedure into a Monte Carlo algorithm. 
The following pseudocode generates random sample paths of the observations Ayi 
and filter pi , sampled faithfully from the probability measure induced by the repeated 
interaction model on the space of observation paths. 
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Time t = lX^ 

Fig. 5.1. (Spontaneous emission) Three typical sample paths of the conditional expectation 
TTi (H) of the atomic energy (top plot) with respect to the integrated photodetection signal Y (l) (bottom 
plot). The atom spontaneously emits exactly one photon at a random time, after which it attains 
its lowest energy and remains there. The initial state chosen was p(x) = Tr[X]/2; in this state 
H = ±?!ti;o/2 have equal probability. The time scale used for the calculation is = 300. 



1. Initialize po- 

2. Z ^ 0. 

3. Repeat 

(a) Calculate 

(b) Sample ^ ^ Uniform[0, 1]. 

(c) If ^ < p+ ^: Aj/i+i ^ +A; Else: A?//+i < A. 

(d) pi+i ^ pi + Api+i{pi, Ayi+i). 

(e) + 

4. Until I = k. 

5.6. Examples. Using the Monte Carlo method developed in the previous sec- 
tion, we can now simulate the observation and filter sample paths for our usual exam- 
ples. In the simulations we have used the initial state p{X) = Tr[X]/2, under which 
the probabilities of the energy attaining its maximal or minimal values are equal. 

Spontaneous emission. In Figure [5TT| photodetection of a spontaneously emit- 
ting atom is simulated. The observation process takes a very simple form: if the atom 
initially attained its maximal energy (which it does with probability ^ under the state 
p), it emits a single photon at a random time. If the atom was already at its lowest 
energy (also with probability ^), the atom never emits a photon. The conditional 
expectation of the energy attains its minimal value immediately after the observation 
of a photon, as at this point we know that the atom has attained its lowest energy. 
Before the observation of a photon, the conditional expectation decays: the longer 
we fail to observe a photon, the higher the (conditional) probability that the atom 
started out with minimal energy to begin with. 

Note that the higher the initial expectation of the energy of the atom, the slower 
the decay of tt; (H) : after all, if the probability of the atom starting out at its minimal 
energy is very small, then we should fail to observe a photon for a very long time 
before concluding that the atom, against all odds, did start off with minimal energy. 
In the extreme case of unit initial probability that the atom has maximal energy 
{p{X) ~ Tr[(T_|-fT_X]), the conditional expectation of the energy is a step function. 
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Fig. 5.2. (Spontaneous emission) Three typical sample paths of the conditional expectation 
TTi{H) of the atomic energy (top plot) with respect to the integrated homodyne photocurrent Y{1) 
(bottom plot). It is difficult to infer much about the atomic energy from the observation process 
using the naked eye, but nonetheless the conditional expectation is revealing (e.g. in the case of the 
blue sample path, the atom initially possessed maximal energy with unit probability) . The initial 
state chosen was p(x) = Tt[X]/2, and the time scale used for the calculation is = 300. 

The latter can be verified directly using the filtering equation for photodetection, 
which shows that A7r;(iJ) = as long as AY{1) = and tti^i{H) = Twoq/2. 

Figure 15.21 shows a simulation of the same system, but now observed using a 
homodyne detection setup. Evidently the way in which information on the atomic 
energy is encoded in the homodyne observations is very different than in the photode- 
tection case; rather than the sudden gain of information when a photon is observed, 
homodyne detection allows us to gradually infer the atomic energy from the noisy 
measurements. The observation process itself is not very revealing to the naked eye, 
but the filter manages to make sense of it. In the case of the blue curve, for example, 
we infer that the atom almost certainly attained its maximal energy value before time 
t ~ 0.5, whereas after time t ^ 1.5 all three paths indicate that most likely the atom 
has attained its lowest energy. 

Dispersive interaction. The dispersive interaction case is quite different than 
the spontaneous emission case. Recall that in this case the energy observable ji{H) = 
H is constant in time. This does not mean, however, that the conditional expectation 
T^i{H) is constant. Whether the energy attains its maximal or minimal value in a 
particular realization determines the mean direction of the phase shift of the outgoing 
light, which can be measured using a homodyne detection setup. As we gain informa- 
tion on the atomic energy, the conditional expectation gradually attracts to the value 
actually taken by the energy H in that realization. 

This behavior is demonstrated in the simulation of Figure 15.31 Each of the filter 
sample paths ni{H) attracts to either +hw[)/2 or —hw[)/2. As the probabilities of 
maximal and minimal atomic energy are equal under the initial state p{X) — Tr[X]/2, 
the sample paths of ni{H) are attracted to ±fta;o/2 with equal probability. It would 
be presumptions to conclude this from the simulation of only three sample paths, 
but if we are willing to believe that ■ni{H) — )■ ±hLUQ/2 with unit probability then the 
result is evident: after all, P{Tri{H)) — P{ji{H)) = p{H) = for all I, so we can 
only have i^i{H) ±hujo/2 if the two possibilities occur with equal probability (the 
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Fig. 5.3. (Dispersive interaction) Three typical sample paths of the conditional expectation 
■7Ti(H) of the atomic energy (top plot) with respect to the integrated homodyne photocurrent Y(r) 
(bottom plot). Though the energy observable is constant in time, the observations contain informa- 
tion on the atomic energy so that the conditional expectation converges to the actual value of the 
energy. The initial state chosen was p{x) = Tr[X]/2, and the time scale used is = 300. 



probabilities change accordingly if we choose a different initial state p) . The fact that 
'!Ti{H) ±fiwo/2 with unit probability can also be rigorously proved, but we will 
postpone this discussion until section [S] 

Note that the behavior of the filter can already be seen by inspecting the observa- 
tion process using the naked eye: though the observation processes are still random, 
the integrated observations have an upward or downward trend depending on the 
value of the atomic energy. This indicates that the actual observation process AF(^) 
is positive (negative) on average if the atomic energy is positive (negative), i.e., posi- 
tive atomic energy leads to an average positive phase shift on the output light, whereas 
negative atomic energy gives rise to a negative average phase shift. 

We have not shown a simulation of the dispersively interacting atom under pho- 
todetection. Though many photons can be observed, a photodetector gives no in- 
formation on the phase of the output light. Hence no information is gained about 
the atomic energy, and the conditional expectation of the atomic energy is constant. 
Evidently the type of detector used makes a big difference in this case. 

6. The reference probability approach. In the previous section we obtained 
explicit nonlinear recursive equations for the conditional expectations 7ri{X). In this 
section we start from scratch and solve the filtering problem in an entirely different 
way. The key idea here is that the Bayes Lemma [33] allows us to express the filtering 
problem in terms of an arbitrary (reference) state; by choosing the "Radon-Nikodym" 
operator V conveniently, we can reduce the calculation to a manageable form. This 
gives rise to a linear recursion for the numerator in the Bayes formula, which we 
denote by ai{X). iriiX) is then calculated as ai{X)/ai{I). 

The reference probability method is widely used in classical filtering theory follow- 
ing the work of Duncan [32] , Mortensen [71] , Zakai [87j , and Kallianpur and Striebel 
[57l[58]. See e.g. [34| for a systematic exposition of reference probabihty methods 
in filtering and control. The corresponding approach in quantum filtering theory, 
adapted below to the discrete setting, was developed in [TO] , 
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6.1. The strategy. Let us begin by outlining what we are going to do. We are 
looking for a way to calculate t^i{X) = P(j;(X)|^). In classical filtering theory it 
has proved advantageous to express the problem, using the Bayes formula, in terms 
of a measure under which the signal (here the atomic observable ji{X)) and the 
observations (here are independent. We will aim to do the same in the quantum 
case. Unlike in classical probability, however, we do not have a suitable quantum 
version of Girsanov's theorem to help us find such a change of measure. We need to 
use a little intuition to obtain the required change-of-state operator V. 

The following simple Lemma will make this task a little easier. 

Lemma 6.1. Let (^,P) be a quantum probability space, <Z £/ a commutative 
subalgebra, and U € be unitary. Define the rotated state Q{X) = ¥{U* XU) on si . 
Then P{U*XU\U*'^U) = U*Q{X\'^)U for any X € "tf'. 

Proof. This is a simple consequence of the definition of conditional expectations. 
Let K e U*'rfU. Then V{U*Q{X\'^)UK) = Q{Q(X \'rf)UKU*) = Q{XUKU*) = 
¥{U*XUK). But as this holds for any K G U*'tfU and as U*Q{X\'^)U S U*'^U, the 
Lemma follows from the definition of the conditional expectation. □ 

How does this help us in our usual setting (^ ® #fe,P)? Recall from section [231 
that % = alg{AZ(i) : i = 1, . . . , (Z = J® A or Z = /(g) (A + A*) for photodetection 
or homodyne detection, respectively), and ^Vi = U{l)*%U{l). For the time being, let 
us fix a time step I and define the state Q\X) = P{U{l)*XU{l)). Then by Lemma 
EH we can write Tri{X) = U{l)*Q^{X\'^i)U{l). Now note that X e ^Sq is of the 
form X (g) J in ^ g) whereas every element of '^i is of the form I ® C . But we 
already know a state under which the initial system and the field are independent: 
this is simply the state P! Hence if we could write Q\X) = f{V{l)*XV{l)) for some 
V{1) G "^/j then we would obtain using the Bayes Lemma 

MX) = u{i) Q {x\%) u{i) = u{iyp{v{irv{i)\%)u{i) ^ ^^'^^ 

Note that we already have by construction Q'{X) = V{U{l)*XU{l)), but U{1) ^ "^Z. 
Hence V{1) = U (l) in the Bayes formula does not work. As we will demonstrate below, 
however, there is a simple trick which we can use to "push" U (l) into the commutant 
'^l without changing the state QK This gives the desired change of state V{1). 

We emphasize that Lemma |6. II is not an essential part of the procedure; we could 
try to find a new reference state S{V*XV) = P{X) and apply the Bayes lemma 
directly to ¥{ji{X)\'3^i) (in fact, the state S can be read off from Eq. (|6.1[) by applying 
Lemma |6. II in the reverse direction). The state P is a very convenient reference state, 
however, as it allows us to use the natural tensor splitting ^ g) Wk and properties 
of the vacuum state (f>'^^ to determine the necessary V{1) with minimal effort. It is 
possible that this procedure could be streamlined in a more general theory for changes 
of state (in the spirit of the techniques widely used in classical probability theory), 
but such a theory is not currently available in quantum probability. 

6.2. Homodyne detection. We first consider the homodyne detection case, 
i.e. Z{1) = A{1) + A*{1). To use Eq. dUI]), we are looking for V{1) € "^Z such that 
Q\X) = F{V{l)*XV{l)). The following trick 46J allows us to find such a V{1) simply 
by modifying the quantum stochastic difference equation for U{1), Eq. (jl.ip . 

Lemma 6.2. Let V be the solution of the following difference equation: 

AV{1) {M+ {AA{1) + AA*{1)) + M° At{l)} ¥{1 ~ 1), ¥(0) = I. 
Then V{1) e and P{V{l)*XV{l)) P{U{l)*XU{l)) for any X e ^ (E) Wk. 
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Proof. The proof relies on the fact that V — p (j)®^ . For simpUcity, let us 
assume that p is a vector state, i.e. p{X) = (u, Xv) for some t; G (we will relax this 
requirement later on). Then P(X) = [v ® <^®^,Xv ® where $ is the vacuum 

vector. The essential property we need is that AA(Z)$®'^ — AA(Z)$®'^ = 0: this 
follows immediately from their definition (as AA(/)$®'' = ^"^^'^ ® a ® <^®^^^ = 0, 
and similarly for AA). Hence using the fact that U{1 — 1) E commutes with 

AA{1) and AA(/), Eq. gT]) gives 

AU{1) V ^ = {M+ AA*{1) + M° At{l)} U{1 - 1) w (g) ^-^^ 

Similarly, any difference equation of the form 

AV{1) ^ {N^ AA{1) + M+ AA*{1) + N- AA{1) + M° At{l)] V{1 - 1) 

satisfies 

AVil) w ® = {M+ AA*{1) + M° At{l)) V{1 -l)v® 

Hence if V{1 - 1) u « = U{1 -l)v® then V{1) w (g) $®*= = U{1) v®^®^. By 
induction Vil) w (g) = U{1) v (g) for any I if V{0) = U{0) = I. Thus 

p{u{iyxu{i)) = (un) V (g $®^ X uii) V (g $®^') 

= {V{1) V (g X Vil) V (g = F{v{iyxvii)), 

regardless of what we choose for and . 

We are now free to choose and N" so that V{1) satisfies the remaining 
requirement V{1) £ "lol. But if we choose = and ~ A/+ as in the statement 
of the Lemma, it follows that V{1) S ^ (g "^(oi C for every /. Indeed, suppose that 
V{1 — 1) S ^ Then V{1) is defined by the recursion as a function of V{1 — 1), 

AZ{1) = AA{1) + AA*{1) G % and M+,M° e which is obviously contained in 
^ g) "^j. The result then follows by induction, and the Lemma is proved. 

It remains to consider the case that p is not a vector state. By linearity we can 
always write p — Tr[y5A'] for some 2x2 density matrix p. But any density matrix 
can be diagonalized (as it is a positive matrix), so that we can write without loss of 
generality p = Xi wiv* + A2 V2V2 where vi is the (normalized) eigenvector of p with 
eigenvalue A2, and V2 is the eigenvector with eigenvalue A2. As the Lemma holds for 
each of the vector states pi.2{X) = (i'i,2, ^^1^2), it must hold for arbitrary p. □ 

The solution of the filtering problem is now remarkably straightforward. 

Definition 6.3. For any atomic observable X € .^0, the unnormalized condi- 
tional expectation 'Ji{X) £ 'S^i is defined as 

ai{X) = U{1)* nV{lYXV{l)\%) Ull). 



Theorem 6.4. The unnormalized conditional expectation (Ji{X) satisfies the 
following linear filtering equation for homodyne detection; 

Aai{X) = (Ji-i{C{X)) At{l) + ai-i{J{X)) AY{1), ao{X) ^ p{X), 

where C{X) is the discrete Lindblad generator and J^{X) was defined in Lemma \5.8i 
Furthermore the noncommutative Kallianpur-Striebel formula holds: 

7:i{X) = ^^ VXG^o. 
cr;(/) 
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Proof. The Kallianpur-Striebel formula is simply Eq. (j6.ip . To obtain the linear 
recursion, we calculate using the discrete Ito rules 

A{V{l)*XV{l)) = V{1 - 1)*£{X)V{1 - 1) At{l) 

+ V{1 - l)*J{X)V{l ~ 1) (AA(0 + AA*{1)). 

Calculating the conditional expectation with respect to , we obtain 

AF{V{l)*XV{l)\%) = ¥{V{1 - 1)*£{X)V{1 - At{l) 

+ FiVil - iyjiX)Vil - {AA{1) + AA*{1)) 

using AA{l) + AA*{l) e But note that % = %-i(^&\g{AA{l) + AA*{1)} and that 
AA{1) + AA*{1) is independent from ^ (g) '^i-i under P. Hence by the independence 
property of the conditional expectation (cf . the last property listed in Table 13. ip and 
the fact that V{1 - l)*XV{l - 1) G (g) for any X £ ^o, we obtain 

AF{V{iyxV{l)\%) ^ F{V(l - l)*C{X)V{l - At{l) 

+ F{V{1 - l)*J{X)V{l - {AAil) + AA*{1)). 

Now multiply from the left by U{1)* and from the right by U{1). Note that 

U{1)* AF{V{l)*XV{l)\%) U{1) = A{U{1)* F{V{l)*XV{l)\%) U{1)), 

because U{l)*CU{l) = U{1 - l)*CU{l - 1) for any C € %-i (see section El]). Fur- 
thermore U{l)*{AA{l) + AA*{l))U{l) = AY{1), and the Theorem follows. □ 

Notice how much simpler the linear recursion is compared to the nonlinear recur- 
sion obtained through the martingale method; for example, this could make digital 
implementation of the linear recursion more straightforward. Nonetheless the two 
methods should give the same answer, as they are both ultimately expressions for the 
same quantity Tri{X). Let us verify this explicitly. Note that 



AiriiX) = A 



MI) 



— - (Ti{X) (ai-i(I) ~cri{l) 



The first term on the right is easily evaluated as 

= 7r,-,{C{X)) At{l) + n-i{J{X)) AY{1), 
whereas we obtain for the second term (using C{I) = 0) 

ai{X) {aiMiy MI)-') = MX) - MX) i^iMJil)) ^Yil). 

Hence we obtain 

A^i{X) = 7r,_i(/:(X)) At{l) + {^,_i(J(X)) - mix) -Ki^^^Ji^I))} AY{V). 

Writing tti{X) — ni-i{X) + Atii{X) and solving for A'Ki{X) gives precisely the non- 
linear recursion obtained in section 15.21 taking into account the identity 

(T^. (T(n)AY(l))-^ I-^iMJ{I))^Y{l) 

ii + nMJii))^Y{i)) =—-^^-—^-^ 

where we have used {AY{1))'^ — At{l). This sheds some light on the seemingly 
complicated structure of the discrete nonlinear filtering equation. 
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6.3. Photodetection. If we naively try to follow the procedure above in the 
photodetection case Z{1) — A{1), we run into problems. Let us see what goes wrong. 
Following the steps in the proof of Lemma 16.21 we reach thepoint where we have to 
choose and N~ so that V{1) € '^1 . But this is impossibleo because we can not get 
rid of the AA* term (as AA*{1)<^'^^ ^ 0), and AA*{1) does not commute with AA(/). 

We are not restricted, however, to using Lemma [6.11 with U = U{1). To deal 
with the photodetection case, suppose that R{1) is some unitary operator of the form 
/ (g) i? in ^ (g) Wk (i.e. it only acts on the field, not on the atom). Define the state 
Q'(X) = ¥{U{l)*Ril)*XR{l)U{l)) and the algebra ^ = R{l)%R{l)* . Suppose that 
there is some V{1) e such that Q'(X) = T{V{l)*XV{l)). Then by LemmaOand 
the Bayes formula, we have for any X e (using R{l)*XR{l) = XR{l)*R{l) ^ X) 

itv\r2>\ uii\TTn\ _ u{i)*R{iyf{v{iyxv{i)\m)R{i)u{i) 



^i{x) u{iyR{iy Q\x\m) R{i)u{i) = 



u{iyR{iy p{v{iyv{i)\^i) r{i)u{i) 



This is precisely as before, except that we have inserted an additional rotation R{1). 
The idea is that if we choose R{1) appropriately, then (AA(/)) i?(Z)* contains 
a AA* term so that we can proceed as in the previous section to find a suitable 
V{1) e M[. A possible choice of R{1) is given in the following Lemma. 

Lemma 6.5. Define R(l) = exp(A(Z) — A*(Z)), I ~ I, . . . ,k. Then R{1) is unitary, 

ARU) = 1^ iAAH) - AA*iiy + H^^^i^ At(/)| R^ 1), 

and we obtain ttie expression 



2 X . 2 - ^ » ' . sin A cos A , . , . . . sin^ A 



R{1) AA(Z) R{iy = (cos^ A-sin^ A) AA(Z)+ ^ (AA(l)+AA*(iy + —- At(l). 

A X'' 



Proof. Unitarity of R{1) is immediate. To obtain the difference equation, note 
that R{1) — exp{AA{l) — AA*{l))R{l — 1), so that we essentially need to calculate 
exp{AA{l) — AA* (l)) = (exp(A(T_ — Acr+));. But it is not difficult to evaluate explicitly 
the matrix exponential of the 2x2 matrix A(cr_ — '^+)'- 



exp(A(o'_ — (T-i-)) = exp 



-A 



/ cos A 
I sin A 



— sin A 
cos A 



cos A / + sin A (<t_ — (T_|_) 



The expression for AR{1) follows directly. To obtain the remaining expression, note 
that R{1 — 1) commutes with AA{1); hence 



R{1) (AA(/)) R{iy = pA^(0-AA-(0 gAA*(0-AA(0 



/ cos A — sin A 
l^sin A cos A 



1 




cos A sin A 
— sin A cos A 



cos^ A sin A cos A 
sin A cos A sin^ A 



from which the result follows immediately. 



^ This is not surprising for the following reason. Applying the spectral theorem to the commuta- 
tive algebra Vi, we obtain a classical measure space (O, JF) on which and P induce different proba- 
bility measures Q and P, respectively. Suppose there exists a V £ such that Q^{X) = ¥{V* XV). 
Then it is not difficult to verify that Q < P with dQ/dP = L{¥{V*V\'ifi)). But note that t(AA(Z)) 
is distributed under Q in the same way as the observation increment AY{1) under P, whereas 
t(AA{i)) = P-a.s. This contradicts Q <C P and hence the existence of V. 
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For brevity, let us call R{1) (AA(/)) R{1)* — AZ{1). We are now in a position to 
repeat Lemma for the photodetection case. 

Lemma 6.6. Let V be the solution to the following difference equation: 

\ sm A cos A cos A / \ A^ 



'^^M+ + cos AAf° + ^^^4^) A<(0| - 1), y(0) 



A ^" ' ' A2 



r/ien y(0 G and ¥{V{1)* XV{1)) = V{U{1)* R{l)*XR{l)U{l)) VAC e » 
Proof Using the quantum Ito rules, it is not difficult to calculate 

A{R{l)U{l)) = |(---) AA(0 + (^cosXM+ - XsmXM° - AA*{1) 



+ i--- )AA{l)+ [ ^A/++cosAAr + ^^5^^^ Atil)^R{l-l)Uil-l). 



We can now follow exactly the procedure in the proof of Lemma 16.21 and using the 
expression for AZ obtained in Lemma [631 the result follows. □ 

We can now obtain the linear filter for photodetection as before. 

Theorem 6.7. The unnormalized conditional expectation 

ai{X) = U{iyR{l)*P{V{l)*XV{l)\m)Ril)U{l), X e ^0, 
satisfies the linear filtering equation for photodetection 

AaiiX) = ai^i{C{X)) At{l) + <n-i{r{X)) (^AY{1) - ^ At{l)^ , ao{X) = p{X), 

where C{X) is the discrete Lindblad generator and T{X) is given by 

1 / A2 



T(X) = I M+*XM+ -X- X^CiX) 

cos^ A \ sin A 

Furthermore the noncommutative Kallianpur-Striebel formula holds: 

(Ji{I) 

Proof. We begin by calculating A{V {I)* XV [l)) using the Ito rules; this is a 
tedious but straightforward calculation, and we will not repeat it here. The result is 

A{V{iyxV{l)) = V{1 - l)*C{X)V{l - 1) At{l) 

+ v{i - irnxwii 1) (^Az(o - ^ At{i)^ , 

where T{X) is given by the expression 

\2 \4 \2 1 

T{X) = M+*XM+ — - M°*XM° — (A^ Af° + Af °*A^) — X. 



sm 



A cos2 A cos^ A cos^ A 
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Using the expression for C{X), the latter is easily transformed into the form given in 
the Theorem. The remainder of the proof is the same as that of Theorem 16.41 □ 
The presence of the sin A and cos A terms in the linear filter for photodetection is 
an artefact of our choice of R{1). In fact, this choice is not unique: many R{1) would 
work and give rise to different linear filters! However, the Kallianpur-Striebel formula 
guarantees that all these linear filters coincide when normalized with the nonlinear 
filter of section [531 To complete our discussion of linear filters, let us show explicitly 
how normalization of the linear filter of Theorem 16.71 gives back, through a series of 
miraculous cancellations, the nonlinear filter for photodetection which we obtained 
earlier through martingale methods. We begin by writing, as in section [6.21 

m-i{I) 

Using the fact that tti{X) — A7r;(X) + 7r(_i(X), this gives explicitly 

A^K^) [I + (Ar(0 - sin^ A)] 

sin^ A 

+ [7r,_i(r(X)) - 7r,^i(X)7r,„i(r(/))] AY {I). 
Using the expression for T(X) and £(/) = 0, we calculate 



= Tn-i{C{X)) Mil) - ^ k/-i(T(X)) - ^/-i(T(/))] At{l) 



+ . 2 5 9 X [ni_,{M+*XM+) - 7r,_i(X) ni_i{M+*M+)] . 
sm A cos^ A 

Hence we obtain 

AMX) [I + ^/-i(r(/)) {AY{1) - sin2 A)] 

= [7ri-i{C{X))-7:i^i{M+*XM+)+7:i^i{X)7:i^i{M+*M+)] At{l) 

cos^ A 

+ [iTi^iiTiX)) - m^iiX) 7ri_i(r(/))] AY{1). 

Next, we claim that 

[/ + ^i_i(r(/))(Ay(0-sin2 A)]"' 

= [/ + cos2A7ri_i(r(/))]"' Ay(/)+ [/-sin^ A 7ri_i(r (/))]"' {I - AY{1)). 

The easiest way to see this is to consider 7r;_i(T(/)) and AY{1) to be classical random 
variables through the spectral theorem; as {AY{1))'^ = AY{1) we conclude that AY{1) 
is a {0, l}-valucd random variable, and the statement follows directly. Using the 
explicit expression for T{X), we find that 

[/-sin^A7r,_i(r(/))]-^ = - 



and that 



[/ + cos2a^,_i(T(/))] ' = -2 



A2 7r^_l(M+*A/^ 
sin^ A 



A2 7r;_i(A/+*M+)' 

Using these expressions, the remainder of the calculation is a straightforward exercise 
and indeed we obtain the expression for At:i{X) as in section [ 
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7. Feedback control. Everything we have done up to this point has been devoid 
of human intervention. An atom sits in free space, emits radiation at its leisure, and 
all we have allowed ourselves to do is to observe the output radiation and to interpret 
it statistically (filtering). In this section we will allow ourselves to manipulate the 
atom in real time; this provides an opportunity for feedback control, which we will 
approach using optimal control theory (section [8|) and Lyapunov methods (section [9]). 

In this section we carefully introduce the concepts required for feedback con- 
trol. As we shall see, some subtleties arise which we explain and address using the 
framework developed in [19], suitably adapted to our discrete setting. Early work on 
quantum feedback control appears in [9l llOj , and a large number of applications have 
been discussed (albeit not in a mathematically rigorous way) in the physics literature, 
see e.g. [30] and the references therein. 

In order to understand feedback, we need to consider how the atomic dynamics 
can be influenced by control actions, and what information is available to determine 
the control actions. We develop these ideas in the following two subsections. In 
subsection 17.11 we describe how a controller can influence the atomic dynamics. This 
will allow us to apply open loop controls, that is, we can apply a deterministic function 
to the control inputs of the system. Though this is not our ultimate goal, it is a helpful 
first step; in subsection 17.21 we will show how to replace this deterministic input by 
some function of the observation history (feedback). This is likely to be advantageous: 
the more information we have, the better we can control! 

Remark 7.1. We take a moment at this point to discuss the usage of the 
term "quantum control" in the literature. Often this term is used to refer to open- 
loop control, rather than feedback control. Such control problems can be reduced to 
deterministic control problems, as we will show in subsection 1 7. II The classical analog 
of this concept would be deterministic control design for the Fokker-Planck equation. 

Our main goal here is to discuss quantum feedback control, where the feedback 
is based on an observations process obtained from the system to be controlled. This 
corresponds to the classical idea of a control system as consisting of a system to 
be controlled (the plant); a sensor which gives rise to an observations process; an 
actuator which allows one to modify the dynamics of the system in real time; and a 
controller, which is a signal processing device that takes the observations as its input 
and produces an actuation signal as its output. 

Sometimes the term quantum feedback control is used in a somewhat broader 
context. One could consider the setup described above in absence of the sensor com- 
ponent. In this case the controller must be treated as being itself a physical system, 
rather than a signal processing device, and the pursuit of this idea leads to a rather 
different theory and applications (see e.g. [84l[86l|52]). This type of feedback is usually 
called coherent feedback or all-optical feedback (in the context of quantum optics), to 
distinguish it from observation-based feedback which we consider here. 

7.1. Open loop control (no feedback). To add a control input to our model, 
recall from section [2751 that the time evolution of an observable X is given by ji{X) = 
U{1)* XU{1), with repeated interaction unitary 

U{1) = lY^=iM{t) = Af (l)Af(2) • • • Af (0, U{0) = /. 
The interaction in every time slice was given by 

M{1) = p-^{3i-iiLi)AA(l)+j,^i{L2)AA'{l)+ji^i{L'2)AA{l)+n^i{L3)At{l)}^ 
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where Li, ^3 G are atomic operators. To add a control input, we simply allow 
ourselves a choice of different M(Z)'s in every time step. To be precise, let us introduce 
a control set il (i.e. these are the values the control can take), and for each m e 11 
we define a set of atomic operators Li{u), L2{u), Lsiu) € S^q. An open loo'p control 
strategy is a fc-tuple u — {ui, . . . , Uk} where Ui e il for all i, and the corresponding 
time evolution — U {I, u)*XU{l, u) is defined by 

Uil, u) = u^) = M^{1, ui)Af "(2, M2) • • • Af ui), UiO, u) = /, 

where the single time slice interaction unitary is given by 

Note that M'^{1, u), for fixed u, depends only on u through ui, . . . , and similarly 
U{1, u) depends only on wi, . . . , uj. We write U{1, u) rather than [/(/, ui, . . . , u;) purely 
for notational convenience; the latter would technically be more appropriate! 
As before, it is convenient to run the definition backwards in time, i.e. 

Uil,u) = l[U,M,{u,) = Mi{ui)Mi^i{ui^i) ■ ■ ■ Miiui), C/(0, u) = /, 
where Mi{u) is given by 

Mi{u) = e-«{iiHAA(i)+i2(«)AA*(0+-L2(«)*AA(0+L3(u)Ai(0}_ (^Yi^i 

These operators are functions of the current control value, and do not depend on the 
full sequence u. The corresponding difference equation is written as 

AU{l,u) = {M^{ui)AA{l) + M+{ui)AA*{l) 

+M-{ui)AA{l)+M°{ui)At{l)}U{l~l,u), C/(0,u)=/. (7.2) 

One could imagine the different controls w G U to correspond to different values of a 
magnetic field which is applied to the atom and can be changed by the experimenter 
in each time step. Another common control input is obtained using a laser whose 
amplitude can be controlled by the experimenter. We will discuss a specific example 
in the context of the discrete model in section [5751 

The question of controller design could already be posed at this deterministic 
level. Though this rules out feedback control (as the observations are a random 
process), open loop controls for quantum systems have already generated important 
applications. For example, optimal control theory has allowed the design of time- 
optimal pulse sequences for nuclear magnetic resonance (NMR) spectroscopy that 
significantly outperform the state-of-the-art for that technology [55] . 

Open loop control design is beyond the scope of this article, so we will not go into 
detail. Let us take a moment, however, to make a connection with this literature. 
Let X g be some atomic operator; then using the quantum Ito rules, we easily 
establish as in section that the expectation of the controlled time evolution j^{X) 
in the case of deterministic u satisfies 

\^^^F{j^_^{C{X,ui))), (7.3) 

where the controlled Lindblad generator is given by 

C{X, u) = M+{uyXM+{u) + M°{u)*XM°{u) + M°{uyX + XM°{u). 
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As P(j"(X)) is linear in X, we can introduce a (deterministic) 2x2 matrix z^" such 
that P(j"(X)) = Tr[i/"X]; then Eq. (|7.3p can be written as a deterministic recursion 
for I/", known as the (controlled, discrete) master equation: 

^ ^ M+iui)i^r-iM+{uir + \^ M°{ui)i^r^,Ar{uir + M°(^^o^/"ii + '^r-iM°{uir. 

The control goal in such a scenario is generally formulated as the desire to choose u 
so that the expectation of a certain atomic observable, or a nonlinear function of such 
expectations, is maximized. But these expectations can be obtained in closed form by 
solving the master equation, so that the control problem reduces to a deterministic 
optimal control problem for the master equation. This is precisely the sort of problem 
that is solved in |59j (in a continuous time context). 

7.2. Closed loop control (feedback). We now turn to the issue of letting 
the control at time step I be a function of the observation history prior to time I. 
Mathematically, it is not entirely clear how to modify the description in the previous 
section to allow for feedback in the repeated interaction model; there are in fact some 
subtleties. It is the goal of this section to clear up this point. For notational simplicity 
we assume from this point onwards that it C M, i.e. that our control input is a real 
scalar (this is not essential; everything we will do can be generalized). 

Let us first consider what we mean by a controller. A controller is a signal 
processing device — a black box — that on every time step takes an observation Ay as 
its input and generates a control signal in il as its output. The most general control 
strategy fj, is thus described by a set of functions, one for each time step: 

^l = {/i, /2(Ayi), /3(Ayi, Ays), • ■ • , /fc(Ayi, . . . , Ay^.i)}. 

Causality is enforced explicitly by making the control in time step I a function only 
of the observations up to that time step. The feedback control /i is called admissible 
if /; takes values in il for every I. We call the set of all admissible controls 

The functions fi encode the input-output behavior of the controller implementing 
the admissible strategy fj, G We now need to hook up the controller to our model 
of the system. As the observation at time i is described by the observable AYi, we 
define the output of the controller at time I as the observable 

u, = /,(Ar(i),...,Ar(/-i))e^^_i. 

Clearly i{ui) is a random variable that takes values in il, precisely as it should be. 

It now remains to close the loop, i.e., to make the time evolution a function of 
the output U; of the controller. Formally, we can proceed exactly as in the open loop 
case; i.e., we define the time evolution jj^{X) = U^{1)* XU^{1) by 

U^{1) = ntiM''(«,u,;) = Af^(l,Ui)Af^(2,U2) • • • A/''(^Uz), t/''(0) = /, 

where the single time slice interaction unitary is given by 

M^^{l,u)=3l^{Mi{u)). 

To make this precise, however, we need to define what we mean by composing the 
unitary operator- valued function M'^{l,u) with the observable u;. We will thus take 
a moment to talk about such compositions. 
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Composition of an operator-valued function and an observable. Let us 

begin with a simple classical analogy. Let u be a random variable in £°°{^l,T,P), 
where J7 is a finite set. Then u takes a finite number of values, and we will suppose 
these lie in a set il. Let M : ii ^ V he a, map from H to some linear space V. We 
would like to define M{u) as a ^-valued random variable. We could do this as follows: 

M(u)(^)= M{u)x{u=u}{oj), (7.4) 

where ranu is the range of u and X{u=ti} is the indicator function on {uj : u(w) = u}. 
This is particularly convenient if we would like to think of f T, P) as being itself 
a linear space; by elementary linear algebra the set of V^-valued random variables is 
isomorphic toV® £°°{i^,!F, P), and the definition of Af (u) above only involves sums 
and tensor products {M{u) (8) X{u=«}) which are naturally defined in this space. 

Though this is obviously not the simplest way of defining composition in the 
classical case, (|7.4|) looks the most natural in the noncommutative context. Let us 
consider the algebra £/ where '^^ is commutative but not necessarily so. We 
can think of jz/ (g) as a linear space of ^-valued random variables; indeed, applying 
the spectral theorem to we find that jz/ «> <^ ~ (g) f °° ( J^) ~ £°°{T; £/). Now 
suppose we are given a map M : ii ^ £/ and an observable u G such that t(u) 
takes values in il. Then it is natural to define the composition Af (u) as an element in 
jz/ (g) in the same way as (|7.4p . This motivates the following definition. 

Definition 7.2. Let £/, ^£ he * -algebras where is commutative. Let il C M, 
M : H and let u G be such that spu = rani(u), the spectrum of u, is a subset 
of il. Then the composition M{\x) £ jz/ (g is defined by 

Af(u)= M{u)P^{u), 
tiesp(u) 

where Pu{u) = {x{l{u)=u}) is the eigenspace projector of u for eigenvalue u. 

In what follows it will be important to understand how the composition A/(u) 
behaves under unitary transformations. Consider jz/®^ C and let J7 be a unitary 
operator in ^. Consider the map Mu : ii U*£/U defined by Mu{u) = U*M{u)U 
for aU M e il. Then U*M{u)U e (g) "^)[/ = C/WC/ (g U*'^U is given by 

U*M{u)U= Y U*M{u)UU*Pu{u)U = 
Mesp(u) 

Y U*M{u)UPu'uu{u) = Mu{U*uU). (7.5) 
ixesp(u) 

Hence unitary rotations preserve the compositions defined above, as long as we re- 
member to rotate both the observable u and the map M{u). 

Controlled quantum flows. Let us return to the controlled time evolution 
Uf'il). Note that Af^(/, •) : U ^ jf_i(-^ g) -^0, whereas € = jf_i('^;_i). 

Hence according to Definition [721 the composition A/(u;) makes sense as an operator 
in the algebra g) «> .-#) C ^i, and by ((73)) 

M{1, ui) = M{1, /(Ar(l), . . . , AY{1 - 1))) = jj'_,{Mi{f{AZ{l), . . . , AZ{1 ~ 1)))). 

For brevity, we will write 



ui = fi{AZ{l),...,AZ{l^l))e%^i 
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SO that M{l,ui) — jj^_j^{Mi{iXi)). In particular, we can now express U'^{1) as 

t^^(0 = l[UiM,{u,) - AfKuOMi-i(ui-i) • • • Mi(ui), C/'^(0) = /, (7.6) 
which gives the controUed quantum stochastic difference equation 

AC/^(0 = {M±(uO AA(0 + M+(u/) AA*{1) 

+M-{ui)AA{l)+M°{iii)At{l)}U''{l-l), ^{0) = !. (7.7) 

The main thing to note is that in order to close the loop in (|7.2[) . the open loop 
controls ui should be replaced by U/ G '^i-i rather than U/ £ "S^i-i- This ensures that 
the corresponding flow jj^{X) depends on the observations history in the right way, by 
virtue of (|7.5p . Other than this subtlety, the treatment of closed loop time evolution 
proceeds much along the same lines as in the absence of feedback. 

The notion of a controlled quantum flow [19j summarizes these ideas in a general 
context. The following definition is a discrete version of this concept and defines a 
generic repeated interaction model with scalar feedback. 

Definition 7.3 (Controlled quantum flow). The quadruple {!d,M,Z,^) s.t. 

1. ilc 

2. Ml : il —^ ^ (g) ./#; , Ml (u) is a unitary operator of the form (j7.1[) V w G il, 

3. Z is an adapted process Z{1) € Wi, I — 1, . . . , A: such that Z{1) is self-adjoint 
and '^i — alg{Z(i) : i = 1, . . . , /} is commutative for every I, 

4- n & ^ is an admissible control strategy, 
defines a controlled quantum flow j['(X) = U {I)* XU ^ {I) . Here U^{1) is given by 
(|7.6p and the corresponding observations process Y^{1) is given by 

AY^il) = Wiiy AZ{1) C/^(/), I = 1, . . . , fc. 

Remark 7.4. For illustrative purposes, we will concentrate in the following 
on the homodyne detection case Z ^ A + A*; the theory for photodetection Z = A 
proceeds along the same lines. The notion of a controlled quantum flow is much 
more general, however, and even allows for feedback to the detection apparatus. For 
example, recall that a homodyne detection setup can measure any of the processes 
e^'^ A + e~'^'^ A* . We could now make (p a function of the past observations (i.e. (p{l) = 
fi{AZ{l), . . . , AZ{1 — 1))), thus feeding back to the homodyne detector; this fits 
within the framework of the controlled quantum flow as it just requires us to use a 
"nonlinear" Z. Feedback to the detector has proven useful for sensitive measurements 
of the phase of an optical pulse, see [531 [3]. We will not consider this further here, 
but everything we will discuss can be adapted to this case as well. We encourage the 
reader to work out the following sections for the feedback scenario of his choice! 

7.3. Filtering in the presence of feedback. Now that we have resolved how 
to model a quantum system with feedback, the next question to be resolved is whether 
filtering still works in this context. Fortunately this is indeed the case, and in fact little 
changes in the proofs. The resulting filters are completely intuitive: one obtains the 
same filter from a controlled quantum flow as one would obtain by first calculating the 
filter with an open loop control, then substituting the feedback law into the filtering 
equation. In this section we will briefly discuss filtering in the presence of feedback 
using the reference probability method. From this point onwards we restrict ourselves 
to the homodyne detection case Z = A + A* . 
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Fix an admissible feedback control strategy ^ ^ K and the corresponding control 
observables U;, U;. We wish to find an expression for 7r['(X) = P(j'['(X)|^^). For this 
to make sense we have to make sure that the self-nondemolition and nondemolition 
properties still hold. If they do not then we did something wrong (recall that these 
properties are essential for a meaningful interpretation of the theory) ; but let us take 
a moment to verify that everything is as it should be. 

Lemma 7.5. The observation algebra '3^^ is commutative (self-nondemolition) 
and £ {'W^'y (nondemolition) for every I — 1, ... ,k and X £ 

Proof. The unitary M/(u/), by construction, commutes with every element of 
'^i-i. Hence U^'{l)* AZ{i) U^'{l) = Ar^(i) for every i < / - 1, and we find that 

[Ar''(i), Ar^(/)] = C/^(0*[A^(i), AZ(0]?7"(0 = 0. 

This establishes self-nondemolition. Nondemolition is established similarly. □ 
To apply the reference probability method, we need a suitable change of state. 

The proof of the following Lemma is omitted as it is identical to that of Lemma 16.21 
Lemma 7.6. Let be the solution of the following difference equation: 

Ay^(/) = {M+{ul){/\A{l) + /\A*{l)) + M°{ul)^^t{l)]V^'{l~l), l/'^(0) = /. 

Then e ^{ and P(T/^(0*Xt/^(0) = f{U^'{l)*XU^'{l)) for any X £ .M ® #'fe. 

From Lemma |3.4[ we immediately obtain the Kallianpur-Striebel formula, 

v) 

(7-8) 

and following the proof of Theorem 16.41 gives the unnormalized controlled filtering 
equation for homodyne detection 

Aaf(X) = <,(/:(X,u,)) A<(0 +af_,(J(X,uO) Ay^(0, = p{X). 

Here we have used the controlled Lindblad generator 

L(X, u) = Af+(M)*XM+(w) + M°{u)*XM°(u) + + XM°(u), 

and we have written 

J{X, u) = XM^{u) + M+(m)*X + }? M° {u)* XM"^ {u) + \^M^{u)*XM°{u). 

As in section we can also normalize this equation; this gives rise to the nonlinear 
controlled filtering equation for homodyne detection: 

A^f(X) = <,(/:(X,uO)At(0 + 

<,(j7(X,u;))-<,(X + A^/:(X,uOX,(J(/,uO) 

{AY'^{l)-nl^_,{J{I,ui))At{l)). 

Finally, we claim that even in the controlled case AY^{1) — Tr^_j^{J'{I,ui)) At{l) is a 
martingale. To show this, it suffices to demonstrate that 7r[^_-^(j7'(J, u;)) Af(/) is the 
predictable part in the Doob decomposition of AF^(Z): but this follows exactly as in 
the proof of the uncontrolled counterpart of this result. Lemma 15.71 
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7.4. The controlled quantum filter. Since 7r['(X) is linear in X, we can find 
a random 2x2 matrix pf such that t(7r['(X)) — Tr[pfX] for every X g In fact, 
the conditional density matrix satisfies the recursion 

pf = ripl„u^,Ayn (7.9) 
where = < = i{ui) - fi{Ay^, Ayi'_,), 

J{p, u) - Tt[J{p, u)] [p + \^C{p, u)) 



r{p,u,Ay) = p + Cip,u) X' + 



l-\^Ti-[J{p,u)]^ 

(Ay-Tr[J(p,u)]A2), 



and 

C{p, u) = M+{u)pM+{u)* + M°{u)pM°{u)* + M°{u)p + pM°{u)\ 
J{p, u) = M+{u)p + pM+{u)* + \^M+{u)pM°{uY + \^M°{u)pM+{u)*. 

To obtain the recursion (17. 9p we can essentially follow the procedure used in section 
15.41 for the uncontrolled case. The only subtlety here is that we need to deal with the 
presence of feedback in terms such as 7r['_-^ (XAf+(ui)) that occur in the recursion for 
TT^^X). The following lemma shows how to do this; it is comparable to the classical 
statement E[/(X,y)|r = y]= E[/(X,y)|r = y]. 

Lemma 7.7. Consider a map X : H — > ^ and an observable ii Cz '^i such that 
u^' = LiWiiyiiUf'il)) takes values in U. Then t(7rf [X(u)]) = Tr[p;X(w''))]. 

Proof. Using Definition 17.21 wc have 

7ri[X{iX)]=7ri ^(")^u(") 

MGsp(u) 

= J2 p[t/^(o*^(«)^u(u)c/^(om] 

«esp(u) 

= J2 P{U''{l)*X{u)U''{l)\m]U''{l)*P^{u)U''{l) 
Mesp(u) 

= ^ ni[X{u)]Pu^(^iyi,u''{i)'{u)- 

M£sp(u) 

Applying l to both sides, the result follows immediately. □ 

It is sometimes more convenient to use the unnormalized form of the filter. As 
(J^ {X) is linear in X, we can proceed exactly as before to find a random 2x2 matrix 
gf such that L{al^{X)) — Tr[gfX] for every X G . The result of Lemma [7J] is 
easily shown to hold also for erf (•)! and we obtain 

= S](^,f_i,<,Ayf), (7.10) 

where 

E(g, u, Ay)^p^ Z(£-, u) + u) Ay. (7.11) 

The conditional density matrix pf can then be calculated as pf — gi['/Tr[pf]. 

Remark 7.8. The filters (|7.9p and (|7.10p define classical controlled Markov pro- 
cesses |61i Chapter 4], where future increments depend explicitly on previous control 
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Fig. 7.1. Cartoon illustrating the operation of a separated feedback control. The only quantity 
stored by the controller at the beginning of time step I is the filter state in the previous time step pi — i. 
The control ui is calculated using a deterministic function gi and is fed back to the atom through an 
actuator. Simultaneously, time slice I of the field interacts with the atom. We subsequently detect 
the field, which gives rise to the observation Ayi. The filter uses Ayi and ui to update the filter 
state to pi . The procedure then repeats in time step I + 1 with pi as basis for the control. 



actions and on a commutative driving process Ay^^. By analogy with the classical 
case, one could even consider the controlled quantum flow as a kind of "quantum con- 
trolled Markov process" , though a precise statement of this concept is not yet used 
in the literature (but see also Remark [43]). 

7.5. Separated strategies. Recall that an admissible strategy /i € ^ is defined 
by a set of feedback functions fJ. = {fi, ■ ■ ■ , fk} that describe how the control values 
depend on the measurement record. Any causal control strategy can be written in this 
form. We already remarked in section [3.41 however, that in many cases the controls 
need not have an arbitrarily complicated dependence on the measurement record — it 
is sufficient to make the feedback at every time step a function of the filter only. In 
other words, many control goals can be attained using only feedback of the form 

uf =gi-i{pl,), l = l,...,k, (7.12) 

where gi takes values in il for any I. Note that such a strategy is in fact admissible as 
p'i^_i is a function of Ay^, . . . , ^y^_i only, so that we could always write gi-iiPi_i) = 
/i(Aj/5', . • . , Ayf_ J for some function 

Both in the case of optimal control (section [5]) and in control design using Lya- 
punov functions (section [5]), we will see that the theory leads very naturally to strate- 
gies of the form l|7.12p . This is highly desirable; not only does such structure sig- 
nificantly simplify the control design procedure, but the resulting controls are also 
much easier to implement in practice. Note that to implement an arbitrary admis- 
sible strategy /z G ^, the controller must have enough internal memory to store the 
entire observation history. One has no other choice, as the feedback signal may be an 
arbitrary function of the control history. On the other hand, to implement a strategy 
of the form (|7.12p the controller only needs to store the current density matrix p'^ at 
any time. The density matrix is then updated recursively using the filter (|7.9p . 

Control strategies of the form ()7.12|) are called separated control .strategies, as the 
controller separates into two parts: a filtering step, which can be solved recursively as 
described above, and a control step which reduces to the simple evaluation of some 
deterministic function of the filter state. The structure of the separated controller 
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is illustrated in Figure 17.11 The ubiquity of this separation structure highlights the 
fundamental importance of filtering in (quantum) stochastic control. 

Let us finally fix the notation. An admissible separated strategy /i is given by 



where : 5 — s- it are functions on the set 5 of 2 x 2 density matrices (self-adjoint, 
nonnegative, unit trace 2x2 matrices). We denote the set of all admissible separated 
strategies by As, and the corresponding feedback signal is given by (|7.12p . 

Unlike in the absence of control (see section [575)1 . the controlled filter pf need not 
be a Markov process as in general fj, € R may have an arbitrary dependence on the 
observation history {p^ is a controlled Markov process in the sense of [HIl Chapter 4]). 
On the other hand, for separated controls p G Ks the process p^ is in fact Markov. 
The following Lemma can be proved as in subsection 15.51 

Lemma 7.9. For separated p £ M.s , the filter pf satisfies the Markov property: 



8. Optimal controL Now that we have described repeated interactions with 
feedback and have obtained filtering equations in this context, the remaining question 
is: how do we choose the control strategy p to achieve a particular control goal? 
In this section we approach this problem using optimal control theory, where the 
control goal is expressed as the desire to minimize a certain cost function. Dynamic 
programming allows us to construct a control strategy and to verify its optimality. 
For an introduction to optimal stochastic control we refer to, e.g., [611 [63l [35l l4l [T4]. 

In optimal control theory, the control objective is expressed in terms of a cost 
function J{p). In this paper we consider mainly the cost function 



This cost is defined in terms of the quantities Q : ii ^ ^ and K E with K 
and Q{u) nonnegative and self-adjoint for all u € H. Large values of these penalty 
quantities correspond to undesirable behavior, while small values correspond to the 
desirable behavior being sought. The penalty Q{u) will be accumulated at each time 
step, contributing to a running cost. For example, it may take the form Q{u) = 
Q + c{u) where Q € ^ penalizes deviation from a desired state, and the function 
c : it — !• [0, oo) penalizes the control effort. K is a terminal cost incurred at the end 
of the time horizon. Ultimately the goal is to find, if possible, an admissible control 
strategy p that minimizes the cost J(/i) 

Remark 8.1. Note that the cost J{p) is defined on a fixed time interval of length 
k. We could also consider this problem on an infinite time horizon, and use cither 
discounted or average cost per unit time extensions of (|8.ip (without the terminal 
cost K). Though the theory of this section extends readily to this scenario, we will 
restrict ourselves to the cost J{p) for concreteness. 

The key step needed to apply dynamic programming is to express the cost function 
in terms of a filtered quantity. This will be discussed further in subsection 18.21 where 
we describe the concept of information state. The following calculation, which makes 



M = {.9o(p),5i(p), • • • ,5fc-i(p)}, 



E^{9ip';)W{p'o, ■ • • , Pf}) - E^(.g(pp|a{pr}) V/ < J < fc. 




(8.1) 
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use of properties of conditional expectations, the spectral theorem, and Lemma 17.71 
expresses the cost (|8.ip in terms of the filter (j7.9p : 



J{y) = P 



.1=1 



^Pbti(Q(uO)l^^-i]Ai(0+P[ji:Wrfe] 

■ 

=1 

k 

^.«JQ(uO])At + .(7rni^]) 
1=1 

k 



.1=1 



(8.2) 



Here is the P'^-expectation on the probability space (ri^,JF^,P^) obtained by 
applying the spectral theorem to "3^^^ . We retain the superscript ^ to emphasize that 
the observation algebra is control-dependent; Y^^{l) and Y'^{1) need not commute! 

The final form of the cost is now a classical cost function for the classical recursion 
(17. 9|) : in other words, we have reduced the control problem to a problem of classical 
optimal control for the classical controlled Markov process (|7.9p . Hence we can apply 
almost directly the ordinary dynamic programming methodology, which we will do 
in the next subsection. We only need to take care to use different probability spaces 
Vl^, for different control strategies, as the corresponding observations need not 
commute; this does not present any additional difficulties, however. 

In order to implement the dynamic programming method, we will use the fact that 
the conditional distribution of Aj/I* given A?/^, . . . , A?/[^j^ can be evaluated explicitly. 
This proceeds in a manner similar to the proof of Lcmma lS.lll but taking into account 
the dependence of the coefficients on the controls using Lemma 17.71 

Lemma 8.2. For any admissible strategy ii € R we have 



whe 



P'^iAyi = ±A I yi'_,] ^ p{Ay = ±A; p^^, <), 



(8.3) 



(8.4) 



depends only on the previous filter state and the control value being applied, and in 
particular is independent of the feedback strategy p,. 

Before we proceed, let us fix some notation. In the following we will encounter 
separated control strategies p,p*,p, € M.s- We will denote the corresponding con- 
trol functions by ^ = {go{p), . . . ,gk-i{p)}, p* = {9q{p), ■ ■ ■ ,gl-i{p)}, and p. = 
{go{p), . . . ,gk-i{p)}- We will also denote by S an arbitrary admissible strat- 
egy, in which case the corresponding feedback process is always denoted by uf . 

8.1. Dynamic Programming. The goal of this section is to give a brief intro- 
duction to dynamic programming. The main results we will need are the dynamic 
programming equation or Bellman equation (|8.6p . which allows us to construct ex- 
plicitly a candidate optimal control strategy in separated form, and the verification 
Lemma 18. 4[ which verifies that the strategy constructed through the dynamic pro- 
gramming equation is indeed optimal. As motivation for the Bellman equation we will 
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first prove the "converse" of the verification lemma: if an optimal separated strategy 
exists, then the Bellman equation follows. The Bellman equation can also be intro- 
duced in a more general context without this assumption, and without any reference 
to separated strategies; this will be discussed briefly at the end of this section. 

Let us begin by considering a separated control strategy £ Rs- A central idea 
in dynamic programming is that one should consider separately the cost incurred over 
different time intervals within the total control interval 0, . . . , fc. To this end, let us 
introduce the cost-to-go 



=1+1 



The quantity Wi{fi, p) is the cost incurred by the control p, over the interval I,. . . ,k, 

given that pf = p. Note in particular that Wo(/i, po) = J{p)- It is straightforward to 
obtain a recursion for Wi{p, p) using the Markov property of the controlled filter: 



Wi{p,p) = E''[Tr[pfQ«+J] At + Wi+,{p, pi^.M = p] 

= Tr[pQ{gi{p))]At+ ^ p{Ay; p,gi{p))Wi+i{p,T{p,gi{p),Ay)), 



(8.5) 



Ay=±\ 



with the terminal condition Wk{p, p) = Tr[pK]. 

The goal of optimal control theory is to find, if possible, a control strategy p* 
that minimizes J{p). For the time being, let us suppose that such a strategy exists 
within the class of separated controls, i.e. that p* G is such that J{p*) < J{p) for 
any p G ^s- It is not unnatural to expect that the same strategy p* also minimizes 
Wi{p,p) for any I and p: if the strategy p* is optimal over the full time interval 
0, . . . , A;, then it should also be optimal over any subinterval I,. . . ,k. If we assume 
that this is indeed the case, then we can simplify the expression for Wi{p*, p): 



Wi{p*,p)= inf Wi{p,p) 



inf 

gi<--;9k 



T^\pQ{gi{p))]At+ J2 Pi^y;P,9i{p))Wi+i{f^,r{p,gi{p),Ay)) 



Ay=±X 



inf 

91 



inf 

91 



inf 

91 



TT[pQ{gi{p))]At+ inf V piAy; p, gi{p))Wi+i{pMP, 9i{p), ^v)) 

91+1, ■■■,9k 



Av=±X 



Tic[pQ{gi{p))]At+ V p{Ay;p,gi{p)) inf Wi+i{p,T{p, gi{p), Ay)) 
^^^^ 



T^[pQ{9i{p))]At+ J2 p{^y;P,9i{p))Wi+i{p*,T{p,gi{p),Ay)) 

Ay=±X 



The following Lemma makes these ideas precise. 

Lemma 8.3. (Dynamic programming equation). Suppose p* e is optimal 
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in i.e. J{n*) < J{^J') for any n G M.s, and define Vi(p) = Wi{p*,p). Then 



Vi{p) 



inf 



TT[pQ{u)]At+ J2 p{Ay;p,u)Vi+i{T{p,u,Ay)) 

Ay=±\ 
Vk{p)^TT[pK], 



< k, 



(8.6) 



for all I and any p that is reachable under p* in the sense that (p^ = p) > 0. 
Moreover, the infimum in (j8.6p is attained at u — gi (p) . 

Proof. In view of the discussion above, it suffices to show that the conditions 
of the Lemma imply that p* minimizes Wi{p,p), i.e., that Wi{p* , p) < Wi{p, p) for 
any p g .^5, I = 0, . . . , and any reachable p. To this end, let us suppose that this 
statement holds for time step / + 1. Then by (jS.Sp 



TiipQ{giip))]At+ J2 p{^y;P,giip))Wi+^{p*,Tip,gi{p),Ay))<Wi{p,p) 

Ay=±\ 

for any p ^ and reachable p. Note that the left hand side is precisely Wi{p' ,p)^ 
where p' is the control strategy that coincides with p at time / and with p* at all 
other times (this follows as Wi+i{p* , p) only depends on the control functions at times 
Z + 1, . . . , /c). We would like to show that the left hand side is bounded from below by 

Wi{p\p)^Tr[pQ{g*,{p))]At+ p{^y;p,gKp))Wi+i{p\np,9Kp),^y)) 

Ay=±X 

for all reachable p. Suppose that this is not the case. We define a new control strategy 
/2 e .^S as follows. At time /, we set 



9i{p) 



gtip) Wi{p*,p)<Wi{p',p), 
gi{p) all other p. 



For any time i ^ l^ we let p coincide with /i*, i.e. gi{p) — g*{p) for i ^ I. Clearly 

Wi{p,p)^TY[pQ(gi{p))]^t+ 

p{Ay;p,gi{p))Wi+i{p*,r{p,gi{p),Ay)) < Wi{p* , p) 

Ay=±\ 

with strict inequality for some reachable p. But then 



Jip* 



E npf-MgLM-i))] + wdp*,pf) 



which contradicts optimality of p*. Hence evidently Wi{p* , p) < Wi{p,p) also for 
time step I. It remains to notice that the statement holds trivially for / = fc as 
Wk{p, p) = Tt:[pK] is independent of p, and hence for any I — 0, . . . ,khy induction. 

Returning to the statement of the Lemma, (j8.6[) can now be obtained as above, 
and the fact that the infimum is attained at g^ [p) follows directly if we compare (|8.6p 
with the expression for Wi{p* , p) given by ()8.5p . Hence the proof is complete. □ 

Let us now consider (|8.6p without assuming that an optimal control p* exists. As 
a backwards in time recursion, this expression makes sense without any reference to 
p* . Vi{p) is uniquely defined for any / = 0, . . . , /c and p: after all, the cost quantities 
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Tt[pQ{u)] and Tr[pii'] are bounded from below, so that the infimum exists and is 
uniquely defined at every step in the recursion starting from the terminal time k. 

Eq. (|8.6p is called the dynamic programming equation or Bellman equa- 
tion. By Lemma I8.3[ we find that if an optimal control fi* e M.s exists, then it 
must be the case that Vi [p) = Wi {fi* , p) and that all the infima in the recursion are 
attained; moreover, in this case an optimal control can be recovered by choosing the 
control functions gi{p) to be minimizers at every step in the Bellman recursion. 

One could now ask whether the converse also holds true. Suppose that in the 
Bellman recursion all the infima turn out to be attained. Then we can define a 
separated control strategy p* by setting 



Qi {p) e argmin 



1:t[pQ{u)]M+ P(A2/;p,u)^i+i(r(p,u,Ay)) 

Ay=±A 



•7) 



The question is, does this imply that p.* is optimal? This is indeed the case, as we will 
show in the following Lemma. Note that this immediately gives a constructive way of 
finding optimal controls, which is precisely the main idea of dynamic programming. 

Lemma 8.4. (Verification and separation). Suppose that all the infima in 
the Bellman equation are attained. Then any control p* defined by (|8.7p is optimal 
in K, i.e. J(p*) < J(p) for any p £ ^. Moreover Vi{p) — Wi{p*,p) for any I and all 
reachable p, and in particular Vo(po) = J(p*)- 

Proof. By substituting (|8.7|) in ()8.6|) and comparing with (18. 5|) . the last statement 
is evident. It remains to show that J{p*) < J{p) for any p G Note that 

Vi{p';)<TT[pfQ{uf^,)]At+ J2 p(Aj;;pr,<+i)^+i(r(pr><+i>Ay)) 

Ay=±X 

by (j8.6p . where we have chosen an arbitrary p (z ^. But by Lemma |8.2[ this is 

mpf)<TT[pfQiul,)]At + E^^m+r{pl,)\yn. 
By recursing the relation backwards from / — fc — 1, we obtain 

k 



Voipo) < 



J2T^T[pl,Q{uf)]At + Tr[p^,K] 



.1=1 



J{p). 



But we had already established that Va{po) = J{p*). Hence the result follows. □ 
Let us reflect for a moment on what we have achieved. We began by showing 
that the cost-to-go for any strategy p* that is optimal in the class of separated 
controls must satisfy the Bellman recursion. Conversely, if the infima in the Bellman 
equation are all attained, we can construct a control strategy p* by solving the Bellman 
recursion {dynamic programming). We then verified in Lemma 18.41 (the verification 
lemma) that the strategy thus constructed is indeed optimal, not only in the class of 
separated controls but in the class R of all admissible controls. Combining these two 
results, we conclude that any strategy that is optimal in .^5 is necessarily optimal 
within the larger class ^ (Lemma 18.41 is also called the separation lemma for this 
reason). This shows that the idea of separated controls is very natural, and indeed 
universal, for this type of control problem. 

Remark 8.5. We have not given conditions under which existence of the infima 
in the Bellman recursion is guaranteed. This can be a delicate issue, see e.g. |15) . 
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However, if il is a finite set the infima are trivially attained. This particularly simple 
case is also the most straightforward to implement on a computer. 

Remark 8.6. The results above are formulated only for reachable p. This is as it 
should be, as for non-reachable p the cost-to-go is defined as a conditional expectation 
on a set of measure zero and is thus not unique. These issues are irrelevant, however, 
to the application of the dynamic programming algorithm and the verification lemma, 
which are the most important results of this section. 

Finally, we briefly remark on the case where the infima in the Bellman equation are 
not attained, ft follows from Lemma l8.3l that there cannot exist an optimal separated 
control in this case. Nonetheless the solution Vi{p) of the Bellman equation, called 
the value function, is still a relevant quantity; it can be shown to be the infimum 
(not minimum!) of the cost-to-go over a suitable class of (not necessarily separated) 
controls. This characterization can be useful in practice, for example, if we wish to 
quantify the performance of sub-optimal controls. 

This approach also provides a different entry point into the theory of dynamic 
programming than the one we have chosen, providing additional insight into the 
structure of the theory. One could begin by proving directly that the value function 
Vi{p), now defined as the infimum of the cost over the time horizon I, . . . ,k, satisfies 
the dynamic programming equation. This does not require us to assume the existence 
of an optimal control, and in particular places no a priori preference on the class of 
separated controls. The only thing that is needed is the fact that the filtering equation 
is a controlled Markov process, which is key to the entire procedure. Indeed, we have 
used this property in an essential way in the form of Lemma [8T2l If the infima in the 
Bellman recursion are attained then we can construct an explicit separated control 
strategy as was demonstrated above; verification then proceeds as we have indicated. 

We will not further detail this approach here, as we do not need these results 
in the following. We refer to [15j for an extensive study of discrete time dynamic 
programming, or to e.g. j35j for the continuous time case. 

8.2. Information States. The key idea that facilitated the dynamic program- 
ming solution to the optimal control problem discussed above was the representation 
of the cost J{p) in terms of the filtered quantity p^*. This is an instance of a general 
methodology involving information state representations |61[ Chapter 6]. An infor- 
mation state is a quantity that is causally computable from information available 
to the controller, i.e. the observation record and previous choices of control values. 
The dynamic evolution of the information state is given by an information state fil- 
ter. Solving a given optimal control problem reduces to expressing the cost function 
exclusively in terms of a suitable information state, then applying the dynamic pro- 
gramming method. The resulting optimal feedback control will be a separated control 
relative to the information state filter. In the previous sections we used the density 
operator as an information state to represent the cost J{p), as indicated in the 
calculation (|8.2p . and the filter for this information state was given by (|7.9p . 

The choice of information state for a given problem is not unique. In the case of 
the cost function J(/i) defined by (|8.I|) . we could also use the unnormalized conditional 
density operator discussed in subsection 17.41 as an information state, with the 
corresponding filter ()7.10p . To see this, we use the reference probability method and 
condition on "^(oi as follows. We begin by defining the state 

and we denote the associated classical state, as obtained through the spectral theorem. 



A DISCRETE INVITATION TO QUANTUM FILTERING AND FEEDBACK CONTROL 63 



by P"^. Under P"^, Ayi, . . . , Ay^ are i.i.d. random variables taking values ±A with 
equal probability, and in particular the law of the process Ay under P"** is independent 
of the feedback control ^ G ^. Then (cf. (|8.2p ) we have 



.1=1 



J{ti) = P XI ^^(' ~ iTQiuiW^'il - 1) Mil) + F^(fc)*i^T/^(fc) 
1=1 

k 

Y,W{1 - l)*Q{ni)V'\l - At{l) + ¥[V''{k)*KV''{k)\'^k 

1=1 

k 

J2 U^'{k)U^'{kyF[V^'{l - l)*g(ui)F^(? - l)|<^,_i]J7^(fc)t/^(fc)* At{l) 
+ U''{k)U''{kyP[V''{k)*KV''{k)\'^k]U''{k)U''{k) 



where the change of state operator V'^{1) was defined in Lemma 17.61 Now recall 
that U^'{k)*V[X\%]U^'{k) = f7^(/)*P[X|<^i]f7^(0, see the proof of the nondemolition 
property in section \Tl5[ Changing to the state P°'', we therefore obtain 



J(/i) = P°'^ 



.1=1 

' k 

.1=1 

' k 



U=l 



where the unnormalized conditional state was defined by 



and is the 



associated density matrix. Using this representation, we can define a value function 
Si{g) and find an optimal control using the alternate dynamic programming equation 



Si{g) = inf 

tiGU 



TT[gQiu)]At+ X p'{Ay)Si+i{^ig,u,Ay)) 

Ay=±\ 

Sk{g) = Tr[£i/v], 



, I < fc. 



where p'^{Ay — ±A) = 0.5. In fact, the optimal control is given by 



h1{g) G argmin 



TY[gQ{u)] At 



Ay=±\ 



p"iAy)Si+i{J:{g,u,Ay)) 



This is a separated feedback control relative to the information state filter (|7.I0p . 

The conditional state (either normalized or unnormalized) is not the correct choice 
for every cost function. In 1981 Whittle fST discovered a different type of informa- 
tion state for optimal control problems with exponential cost functions (risk-sensitive 
control), which are not solvable using the standard conditional states. Instead, the 
filter for the corresponding information state depends explicitly on quantities defining 
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the cost function, and the optimal feedback control is separated relative to this filter. 
Following [S3] we now explain this briefly in the context of this paper. 
In the quantum setting, a risk-sensitive cost can be defined by 



C{k) 



C{k) 



where AC(/) = f jf_i(g(u;_i)) C(/ - 1) At(0, C(0) = / defines the "exponential" 
running cost and 6* > is a fixed real parameter (the risk parameter), 
define [/^(Z) — U'^{l)C{l) (which is not unitary in general), so that 



Let us now 



We can now proceed as in the previous part of this section to express the control cost 
in terms of an unnormalized filtered quantity. The corresponding filter is not obtained 
from the usual unitary U^{1), however, but from the modified operator U^{1). We 
can obtain a change-of-state operator V^{1) as in Lemma [7.61 which gives rise to an 
information state filter that depends explicitly on the running cost Q{u). We can 
subsequently express J^(/i) in terms of this filter, and the optimal control problem 
can then be solved using dynamic programming. We leave the details as an exercise. 

In classical stochastic control, the risk-sensitive cost is known to possess improved 
robustness properties compared to the usual cost J(/i); in particular, as the risk 
parameter 9 increases, the optimal performance becomes less sensitive to the details 
of the underlying model. To what extent these advantages carry over to the quantum 
case remains to be explored (but see [5T]). 

8.3. Example. We will give a numerical example of dynamic programming for 
a particularly simple system — a controlled version of the dispersive interaction model. 

The controlled quantum flovif. We consider again the dispersive interaction 
model of section [221 but now we add a control input. The controlled repeated inter- 
action matrices are now given by 



Li{u) = 0, L2iu) = iaz 



Lsiu) = iu(a^ ~ a J). 



Such an interaction can be realized in certain systems by applying a magnetic field 
of strength u to the atom, see e.g. Figure II. 1[ in addition to the usual dispersive 
interaction with the electromagnetic field. In principle any m g R is admissible, i.e. 
we should take il = R. As we will be evaluating the dynamic programming recursion 
numerically, however, it is more convenient to discretize the admissible controls, i.e. 
we choose il = {— i^e, — (K — l)e, . . . , Ke\ C K where e is the discretization step size 
and 2K + 1 is the total number of admissible controls. 

Our first task is to evaluate the matrices M^'+'~'°(u) in the controlled difference 
equation. Let us calculate explicitly 



Mi(u) = exp(cr^(AA*(0 - AA(0) + u(g+ - a^)At{l)). 



Writing 



B{u) = (g) A (cr+ - a-_) + u(cr+ - cr_) (g> X^I 








A 


M 










-A 


-A 








A^M 





A 


-A^u 


/ 
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we calculate the matrix exponential 



/ cos{Xw{u)) uX ^'"(^"'(")) 



\ sin(Xw(u)) 
iA 7 — \ 

w{u) 

sin(Aiu(u)) 
w{u) 





v(u) 

cos(Aw(u)) 


w{u) 



sin{Xw(u)) 
w(u) 



cos(Au'(u)) 

sin{\w{u)) 



\ 

sin(\w (u)) 

w{u) 
sin{Xw(u)) 
w(u) 

cos{\w{u)) J 



where we have written w(u) = + Hence we obtain 



, sm( \w(u)) ,. , ,,,,, 

^'(") = ^ , \ (AA*(/) - AA{1)) 



Xw{u) 



cos(Aui(u)) — 1 sin(Aw(M)) 



A2 



Xw(u) 



At{l)+I. 



We can now immediately read off the coefficients in the quantum stochastic difference 
equation for the controlled dispersive interaction model: 



An invariant set. The solution of the filtering recursion is always a density 
matrix: recall that t(7r['(X)) = Tr[p['X] is a positive map with 7r['(/) = 1, so that Pi 
must be a positive matrix and have unit trace. The the goal of dynamic programming 
is to calculate the feedback function for any time step / and density matrix p. 

Unfortunately, even the space of 2 x 2 density matrices is rather large; discretization 
of this space, as one would need for computer implementation of the dynamic pro- 
gramming recursion, would require a tremendous number of discretization points. For 
this reason dynamic programming is computationally expensive, prohibitively so in 
moderate- to high-dimensional systems where optimal control theory often plays the 
role of a benchmark rather than a practical solution. In such cases, one is compelled 
to settle for control designs that are sub-optimal, i.e. they do not minimize a cost 
function. We will briefiy explore one such approach in section [S) 

The simple example treated in this section, however, has a feature which signif- 
icantly simplifies the implementation of dynamic programming. As we will demon- 
strate shortly, there is an invariant set of density matrices which is parametrized by 
points on a circle: i.e., if we start the filter somewhere on this circle, it will always 
remain there. This reduces the dynamic programming algorithm to the calculation of 
the feedback function gi{p) on the circle. With a sufficiently accurate discretization 
of the circle, this problem can be solved numerically without too many complications. 
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To find the desired invariant set, note that 



„„, , sin(2Au'(u)) , , sin^fAwfw)) 
C^ia,, u) = \ , ; u a+ + a_ 



AW(U) 

„„, X sin(2Afi;(w)) 2sin^(Aw(u)) , , ,„ „, 

C>'{a+ + cj.,u) = i-^-A-^wa, ^ " {a++u-), (8.8) 

, X 2Msin^(Aw(M)) 

w{uY 

tp.it „.x _ sin(2Aw(u)) _ , 2m sin^(Au;(M)) 
Xw{u) w{u)^ 

Hence evidently the recursion for 7:^{a.^) and 7r['((T+ + (t_) forms a closed set of 
equations. We claim furthermore that TT^{azY + + <^-)^ = ^ fo^' ^■ll ^) if this is 

true for I = 0. The algebra is a little easier if we consider the unnormalized version; 
using the discrete Ito rule, we calculate 

A(af(X)2) = 

{2ar_,(X)ar_,(£''(X,<)) + AV-i(-C''(X,<))2+<_,(J''(X,<))2}At(0 

A tedious but straightforward calculation shows that 

A[ar(a,)2+ar(a++a_)2-ar(/)2] = [ar_,(a,)2+ar_,(a++a_)2-ar_,(/)2] At{l) 

where 'd^.x is a complicated function of u and A. Hence if ttq (cr;2:)^+7rQ ((t++<t_)^ = I, 
it follows that 7rf (ct^)^ + 7rf (cr+ + o--)^ = / for aU L 

Now let p be any density matrix. The remaining insight we need is that if we are 
given x,z such that Tr[p(tT+ + a-)] = x, TT[paz] = z, and + = 1, then this 
uniquely determines p. To see this, let us write without loss of generality a: = sin^ 
and z = cos6. Using the constraints Tr[p{a+ + ct-)] = sin 6*, Tr [pcTz] = cos 6*, Trp = 1 
and p = p*, we easily find that 

/ i + i cos 6* i sin ^ + « ^ « 

P=\? ■ a -a \ 1 a) for some /? € M. 
\^s>m.d — ip | — ^cos^y 

But we can explicitly calculate the eigenvalues of this matrix as ^(1 ± (1 + 4/3^)^/^), 
so that the remaining requirement for the density matrix p > implies that /3 = 0. 
Hence we conclude that the "circle" of density matrices 

parametrized by the angle is invariant under the filtering equation for our controlled 
quantum flow in the sense that G for all / if po S 5^. We can thus restrict the 
dynamic programming recursion to this set, which yields a feedback control law on 
the circle of the form gi(&). 
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Fig. 8.1. The optimal control function g1{0), for selected I, obtained by applying the dynamic 
programming algorithm to the cost function Eq. I8.9il . During the major portion of the control 
interval, the feedback function is essentially equal to the first plot; the control passes through zero at 
the state of maximal energy, and has a break at the point of minimal energy. Close to the terminal 
time, the optimal control varies a little to accommodate the terminal cost. 



Dynamic programming. We are now in a position to solve the dynamic pro- 
gramming algorithm numerically. To this end, we have discretized the circle into a set 
of 10^ equidistant points, and we have discretized the control set il into 400 equidis- 
tant points in the interval [—10, 10]. As in the previous simulations, we have chosen 

= 300 and a terminal time of 3 (i.e. k = 900). 

As a first control goal, suppose we would like to maximize the expected energy, 
i.e. we would like to drive TilpfcTz] to -1-1. To this end, we use the cost (|8.ip with 

Q{u)^Cu^+D{I~a,), 

K = I-a, (8.9) 

The first term in Q{u) penalizes the use of large control strengths, which is necessary 
in any practical feedback loop. The second term tries to minimize 1 — Tr[p['crz] during 
the running time of the system, whereas the terminal cost tries to minimize the 
terminal value of 1 — Tr[pfaz]- The constants C,D > determine the relative weights 
attributed to these control goals. As an example, we have chosen C = 0.25, D — 5. 
The corresponding optimal feedback function gi{0) is plotted for several times / in 
Figure 18.11 Note that the control passes through zero at the state of maximal energy 
9 = 0, whereas the break at 9 = tt drives the system away from the state of minimal 
energy. This is not unexpected, as both 9 = {0,7r} are fixed points of the filter (see 
section [5^ : the break in the control destabilizes the undesired minimal energy state. 

As a second example, suppose we would like to drive the system close to the point 
where Trfp^'uz] = Tr[/3['cr2:] = 2~^/^, i.e. 9 = 7r/4. This point is uniquely characterized 
by Tr[p['((72 -|- a^)] = V^, so we choose 

Q{u) = Cu^ + D{I - X), 

K = I-X (8.10) 

where X = 2^'^/'^{(Tx + Uz) and, for example, C = 0.25 and D = 5. The optimal 
feedback function for this case is plotted in Figure 18.21 Once again, the feedback 
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Fig. 8.2. The optimal control function /or the cost Eq. \8.10\l . As before, the feedback 

function is essentially equal to the first plot durinq the major portion of the control interval, and 
the control passes through zero at the target state. In this case, however, the control function is not 
discontinuous at the point of the circle opposite to the target state. 




Fig. 8.3. The cost J(/x*), as a function of po £ , for the optimal control strategy p," min- 
imizing the cost (a) Eq. S8.91) . and (b) Eq. f8.10\} . The cost is simply obtained as the final value 
function in the dynamic programming recursion, J{fi*){po) = Vo[po]- case (a) the cost is zero at 
9 = 0: after all, the target point is a fixed point of the filter, so no cost is accumulated if we start at 
the target point. This is not the case in (b) where the target point is not a fixed point. 



function passes through zero at the target point 9 = 7r/4. However, note that the 
function is no longer singular; as the opposite point on the circle is not a fixed point 
of the filter in this case, it is evidently more efficient to let the filter drive itself toward 
the target point without expending additional control effort. 

Finally we have plotted the optimal cost J(/i*), as a function of the initial state p G 
S^, in Figure lSTSl The cost is easily obtained as a byproduct of dynamic programming, 
as it is simply given by the value function at time zero. Note that in the case of our 
first example, zero cost is possible: as the target = is a fixed point of the filter, 
no cost is accumulated if the system is initially in its state of maximal energy. In our 
second example, on the other hand, this is not the case: even if we start at the target 
state 9 = 7r/4, the filter will fluctuate around this point and a total cost (just under 
J[/x*] = 2) is accumulated over the control interval. 
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9. Lyapunov control. Many control design methods other than optimal control 
have been developed and are applied in control engineering, for a variety of reasons. 
Indeed, the applicability of dynamic programming is limited by computational com- 
plexity for use in low dimensional problems, or in situations where explicit solutions 
are available. Among the alternative methods are Lyapunov design methods, and the 
purpose of this section is to demonstrate their use in our quantum context. 

Consider for example the following scenario. Recall that the dispersively interact- 
ing atom, when untouched by the experimenter, has constant but unknown energy; 
our control goal is to drive the energy to a particular value of our choice, say +hujQ/2. 
In the absence of sharp time constraints this problem is not uniquely expressed as 
an optimal control problem: what cost should one choose? All we want to achieve 
is that V{ji{H)) hujQ/2 for large I, a much more modest control goal than the 
absolute minimization of a particular cost function. Hence dynamic programming is 
unnecessarily complicated for the problem at hand, and we can resort to a simpler 
method for control design. In this section we will design a controller for the control 
problem described above through a very simple Lyapunov function method, similar 
in spirit to 79J. The simplicity of the method is its chief feature: not only will we 
find a feedback function that works, but the resulting feedback function is also of a 
very simple form and is easily implementable. 

As before, we can easily express the control goal in terms of the filter: P{ji{H)) = 
P(7r;(iJ)) huJo/2 for large I. We will design a separated control that achieves 
this goal. The Markov property of the filter with separated controls is key to the 
procedure: it allows us to treat the problem entirely at the level of the filter, without 
reference to the original repeated interaction model. 

9.1. Lyapunov stability. The main tool is the following Lyapunov theorem for 
Markov processes, taken directly from [63| . 

Theorem 9.1 (Lyapunov theorem). Let xi he a (classical) Markov process. 
Suppose there exists a nonnegative (Lyapunov) function V{x) > that satisfies 

E{AV{xi)\a{xi-i}) = E{V{xi)\a{xi^i}) - V{xi^i) = -k{xi^i) < V/, 

where k{x) > is another nonnegative function. Then fc(x;) — > as I oo a.s. 
Proof. Fix some n > 0. As V{x) is a nonnegative function, clearly 

V{xo) > V{xo)--E{V{xr.)). 



But by the condition of the Theorem, we obtain 



E{V{x„)) - E(F(x„_i) - fc(a;„_i)) = E(F(x„_2) - Hxn-2) - K^n-i)) = ■■■ 
Iterating this procedure, we obtain 

As this holds for any n and as l^(xo) < cxd, the right hand side is finite for any n. But 
PiHxi) > e) < - E ^ k{xi) < 00 Ve > 

l=Q ^ \l=Q J 
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by Chebyshev's inequality, so the Borel-Cantelli Lemma gives k{xi) a.s. □ 

Remark 9.2. Previously we did everything on a fixed time horizon I — 0^ . . . , k. 
Now, however, we are considering what happens as ^ ^ oo, so technically a little more 
care is needed. The reader should convince himself that not much changes in this case. 
In particular, there is no need to deal directly with an infinite-dimensional algebra 
Wcx, = corresponding to an infinite number of time slices of the electromagnetic 

field. Rather, one can use the natural embedding i : ^ ®Wk^ ^ ® ^i,^) = 

X ® I and the sequence of probability spaces ^ C obtained by applying the 

spectral theorem to the corresponding observation algebras, to give meaning to infinite 
time limits without functional analytic complications (provided the time step A and 
repeated interaction model is fixed). We leave the details as an exercise. 

Before we design the controller, let us verify that without control it is indeed the 
case that niiff) — > -^hio^^j^ a.s. This already demonstrates the Lyapunov theorem. 

Lemma 9.3. //itf — 0, then Tr[pf(T2] ±1 with unit probability. 

Proof. We will use the Lyapunov function V{pf) = 1 — Tilpfaz]"^ , which is 
nonnegative and is zero precisely when Tr[/9[' cr^] = ±1. Using the filtering equation 
with = and Eq. (|8^ . we find that 



sin(2A) V{pl,) 

ATr[p^a,] = . 2,r,y,,r^ ^ n ^^F^ 

A cos^(2A) + sm (2A)\/(pJl^) 

where is the innovations process. Using the quantum Ito rule 

sin^(2A) Vjpti? 

^^^^^^ ^ - A^ (cos^(2A) + sin^(2A)l^(pr„,))^ ^^"'^ +i---)^y^- 

But E''(Ay['|^^^-^) = by the martingale property, and furthermore 

E^((Ayr)2|^^1 J = (1 - A^ TT[J{pl,r) At = (cos2(2A) + sm\2X)Vipl,)) At, 

where we have used iAY{l)f = (•••) AY{1) + (AF(/))^ - i^C{l)f as in the proof of 
Lemma \5M Hence we find 



mAvipf)W{pl,}) = - ""'f^ .(,,Z^''€xW( ^ ^ 

A^ cos^(2A) + sm (2A)\/ 
The Lemma now follows from the Lyapunov theorem. □ 

9.2. Construction of a Lyapunov controL We now turn to the control de- 
sign. We wish to find a (time-invariant) feedback control uf = g{p'^_i) so that 
P{TT^{az)) +1 as I oo. The way we approach this problem is to use a trial 
Lyapunov function V{p) without fixing the control g. By inspecting the expression 
for {AV {p'j^)\a{p'^_j^}) , we can subsequently try to choose g so that this expression 
is nonpositive and is zero only at (or in a small neighborhood of) the point p G 
where Tr[/5cr^] = +1. The desired result follows by dominated convergence. 

Let us implement this procedure. Choose the trial Lyapunov function V{p) = 
1 — Tt[p(Tz]. Using the filtering equation with = g{pf_^), Eq. (|8.8p . and the 
martingale property of the innovations, we find that 



W{AV{p\')\a{pl,}) = 
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As a first attempt, consider the following feedback function: 

r -1 ifTr[p(a+ + a_)] <0, 

■'^P>-\ +1 if Tr[p(a+ + a_)] >0. 

Then 

E>^{AV{pf)\a{pl,}) sin(2Ml)) , ^ , 2 sin^(Au;(l)) ^ ^ ^ 

Keeping in mind that pf E for all I, clearly there exists some S > such that this 
expression is strictly negative for Tr[p['_ j^cr^] < 1 — S; in fact, 5 is the solution of 

sin(2A«;(l)) , 2 sin^(A7^(l)) 
MT^^('-^)'^+ ^(1)2 (l-5) = 0. 

To keep E^(Ay(pf )|(t{pJ'_-^}) nonpositive everywhere, we now modify the control 
function to turn off the feedback in the set Tr[/9f_ -^ctz] > 1 — 6: 

( if Tr[p(7,] > 1 - (5, 

9ip)^{ -1 if Tr[p(CT+ < and Tr[pCT^] < 1 -(5, (9.1) 

[ +1 if Tt[p{<7+ + aJ)] > and TrlpaJ < 1 - S. 

This ensures that E''(A1/(/9[')|(t{p['_^}) < 0, where the equality holds only when 
Tr[p['_j^(Tz] > 1 — 5. The Lyapunov theorem then guarantees that Tr[pf(Tz] converges 
to the set [1 — 6,1]: i.e. liminf;_+oo Tr[p['cr^] > 1 ~ 6 with unit probability. For 
A <C 1, the threshold 6 will be very small and hence the Lyapunov theorem guarantees 
convergence to a tiny neighborhood of the desired control goal. For example, with 
the choice A~^ — 300 which we have used in the simulations, 6 w 5.6 x 10^^. 

Remark 9.4. We can do better than proving convergence to a small neigh- 
borhood of the target point. In addition to the convergence k{xt) 0, the invariant 
set theorems [53] tell us that Xt itself converges to the largest invariant set contained 
in {x : k{x) = 0} (the convergence is in probability, but this is easily strengthened 
to a.s. convergence). As the only invariant set inside the set Tr[/5fTz] > 1 — 5 is the 
target point Tr[yOcr2] = 1, convergence is guaranteed. A full discussion of the required 
theorems is beyond the scope of this article, and we refer to [63^ for further details. 

To illustrate the effectiveness of Eq. (|9.1I) , we have used the Monte Carlo method 
to simulate in Figure [9?T] several sample paths of the controlled filter. It is immediately 
evident that the control goal was attained for the five sample paths plotted, and the 
Lyapunov theorem guarantees that this is indeed the case for any sample path. 

The discussion above is mainly intended as an illustration of an alternative to 
optimal control theory, and we have only used the Lyapunov theory in its simplest 
form. Lyapunov methods can be extended to a large class of systems, see e.g. |70) . 
whose treatment using optimal control would be out of the question due to the high 
dimensionality of the state space. Some control goals other than strict convergence 
can be treated using similar methods; for example, we could try to find a control law 
so that a quantity such as E(|Tr[p['crz] — a\), for some a £ (—1, 1), becomes small 
(but usually nonzero) as Z ^ oo, see e.g. |5S]. Beside the restricted set of control goals 
that can be treated by Lyapunov methods, a drawback of such methods is that there 
is no general recipe which one can follow for the design of controls; the choice of a 
proper control and Lyapunov function has to be investigated on a case-by-case basis. 
Ultimately the control goal of interest and the available resources determine which 
method of control design — be it optimal control, Lyapunov methods, or some other 
approach — is most suitable for the problem at hand. 
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1.5 

Time t = lX^ 



Fig. 9.1. Five typical sample paths of the conditional expectation tti{H) of the atomic energy 
(top plot) for the controlled quantum flow on the circle with separated feedback law Eq. 19. Note 
that iTi{H) = +/ja;o/2 is always attained. The bottom plot shows the feedback signal = g{Pi_-i). 
The initial state was p{x) = (3>, X<E>), p G , the time scale is X~'^ = 300, and 5 = 6 x 10~^. 



10. References for further reading. We have come to the end of our expo- 
sition on the filtering and control of discrete quantum models, though in many ways 
we have only scratched the surface of the theory of quantum probability, filtering and 
stochastic control. The goal of this final section is to provide some entry points into 
the literature for further reading. 

Quantum probability. We have used quantum probability theory in its sim- 
plest finite-dimensional form, where only finite state random variables are available. 
In practice such a description is very restrictive, and one needs a theory that ad- 
mits continuous random variables. The theory of operator algebras (e.g. Kadison and 
Ringrose [Ml [55] or Bratteli and Robinson [21]), particularly Von Neumann algebras, 
provides the proper generalization of the concept of a quantum probability spaces to 
the infinite-dimensional context. For an accessible introduction to quantum probabil- 
ity we refer to the excellent lecture notes by Maassen [551 . See also (THl [1^ and the 
references therein. The required background on functional analysis can be found e.g. 
in the textbook by Reed and Simon [74] . 

Readers with a background in physics might wonder why we did not use the 
projection postulate in this paper. We do not need this postulate to develop the 
theory, and in fact the projection postulate can be obtained as a special case within 
the quantum probability framework used in this paper. See [201 for an example. 

Quantum noise and stochastic calculus. The three discrete noises — we have 
suggestively denoted them as A{1), A*{1) and A{1) — are replaced in the continuous 
time theory by the standard noises At (the annihilation process), A^ (the creation 
process) and At (the gauge process). For each time t, At, At, At are defined as opera- 
tors acting on a Hilbert space which is known as the Fock space, and the latter plays 
a central role in the theory of quantum noise. The processes [At -\- Al)t>o or {At)t>o 
each generate a commutative algebra; the spectral theorem allows us to interpret the 
former as a diffusion process (a Wiener process in the vacuum state), and the latter 
as a counting process (a Poisson process in a so-called coherent state, and a.s. zero in 
the vacuum). A brief introduction is given in e.g. [20] . 
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Based on these processes, one can now proceed to define quantum stochastic inte- 
grals and obtain the quantum Ito rule: dAf dA^ = dt, dAt dAl = dA^, dAt dAt = dAt, 
and all other combinations are zero. Next one defines quantum stochastic differen- 
tial equations, which provide "quantum noisy" versions of the famous Schrodinger 
equation and form good models of actual physical systems, particularly those used in 
quantum optics (see [36] for a physicist's perspective). An accessible starting point 
in the literature on quantum stochastic calculus are Hudson's lecture notes [49]; see 
also the original article by Hudson and Parthasarathy [SO] and Parthasarathy's book 
frl\ . Several other approaches to the theory have been developed since; some of these 
appear in the book by Meyer [69] and the lecture notes of Biane [16] . 

Physical models. The physical theory that describes the interaction of light 
(the electromagnetic field) with matter (atoms, molecules, etc.) is called quantum 
electrodynamics; see the book by Cohen- Tannoudji et al. [25j . Very often, and in 
particular for the purposes of filtering and control, quantum electrodynamics does 
not directly provide a usable model; these models are non-Markov. In many systems, 
however, the Markov approximation is quite good, and applying this approximation 
to models from quantum electrodynamics gives rise precisely to quantum stochastic 
differential equations. The Markov approximation can be pursued at different levels of 
rigor, ranging from the ad-hoc "whitening" of the noise which is common in the physics 
literature [35] to rigorous Wong-Zakai type limits [331 [I] • As in classical probability 
theory it is also common to begin with a Markov model in the form of a quantum 
stochastic differential equation, which is justified either on phenomenological grounds 
or through a separate modelling effort. 

The modelling of the detectors, e.g. photodetection or homodyne detection, is 
another important topic. The operating principles of optical detectors are described 
in many textbooks on quantum optics, see e.g. [501 [SZ]- In the quantum stochastic 
context we refer to e.g. [36', or to the more mathematical perspective of Barchielli [7]. 

Continuous limits of discrete models. The discrete models which we have 
used throughout the paper are conceptually very similar to the more realistic contin- 
uous models. There is more than a conceptual resemblance, however: these models 
in fact converge to corresponding continuous models as we let the time step — > 0. 
Though there is no particular physical intuition behind this convergence except the 
one we have already given — that small time slices of the field approximately contain 
at most one photon — the convergence demonstrates that even results obtained using 
these highly simplified models have some relevance to the physical world. 

A simple example of this convergence is easily verified. Recall that the classical 
stochastic process x{l) = i.{A{l) + A*{1)), with the measure induced by the vacuum 
state, defines a symmetric random walk with step size A. First, let us express this 
process in terms of time t = IX^ rather than the time step number I, i.e., Xt = 
x([t/A^J). Taking the limit (in law) as A ^ gives, by the usual functional central 
limit argument (e.g. [40l pp. 452-]), a Wiener process. But t{At + A^) is a Wiener 
process, so we see that A{1) + A*{1) converges to At + A^ at least in a weak sense. 

Similarly, most of the expressions given in this article converge as A —> to 
their continuous time counterparts. In particular the difference equation Eq. (j4.ip . 
the filtering equations, etc. converge to (quantum) stochastic differential equations. 
Even the discrete Ito table "converges" to the quantum Ito rules dAt dA^ — dt etc., 
if we formally set A A — > dA, A A — > dA, AA* dA*, At dt, and terms such as 
X^AA (this was noticed in [S]). Evidently the discrete theory mirrors its contin- 
uous counterpart quite faithfully, which is precisely the motivation for this article. 
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Convergence of the discrete models is investigated by Lindsay and Parthasarathy 
[5i] . Attal and Pautrat and Gough [55]. These articles demonstrate convergence 
(of various types) of the solution of the difference equation Eq. (|4.ip to the solution of 
a corresponding quantum stochastic differential equation. The convergence of discrete 
filters to their continuous counterparts is investigated in Gough and Sobolev [44] . 

Conditional expectations and quantum filtering. The concept of condi- 
tional expectations in Von Neumann algebras is an old one, see e.g. [78 . However, a 
much more pragmatic notion of conditional expectations on which the definition we 
have given is based first appeared in the pioneering work of Belavkin, see e.g. [12] . 
The necessity of the nondemolition property, which is the key to the development 
of both filtering and feedback control, appeared already in [5] . This opens the door 
for the development of nonlinear filtering, which is done in detail in the difhcult and 
technical paper [12, using martingale methods. The reference probability approach 
to quantum filtering, using the Bayes formula, appears in [T^j (but see also [TT[ BT] ). 
An introduction to quantum filtering is given in [20] using both methods. 

Feedback control. An early investigation on optimal feedback control in dis- 
crete time can be found in Belavkin's 1983 paper [9]. The Bellman equation for 
continuous optimal control appears in [10^. A recent exposition on quantum opti- 
mal control, and in particular quantum LQG control, is given in |33] . see also |17] . 
and see [311 [30] for a physicist's perspective. The separation theorem for the con- 
tinuous time case appears in [19] . A different type of optimal control problem, the 
risk-sensitive control problem with exponential cost, is investigated in [53] . Finally, 
Lyapunov function methods for control design are investigated in [79] 170] . 

Many applications of quantum feedback control have been considered in the 
physics literature; we refer to the introduction for references. Note that the "quantum 
trajectory equations" or "stochastic master equations" used in the physics literature 
are precisely filtering equations, written in terms of the innovations process as the 
driving noise rather than being driven by the observations. The notion of filtering 
has only recently gained acceptance in the physics literature, though the relevant 
equations had already been developed by physicists. The reader should beware of the 
discrepancy in terminology and interpretation between the physical and mathematical 
literature. However, as the Markov property of the filter allows us to pretend that 
the filter itself is the system to be controlled, most of the control schemes found in 
the physics literature make sense in both contexts. 
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